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REMARKS 

Claims 234-247 and 466-479 are all the claims pending in the present application. 
Claims 1-233, 248-465 and 480-509 have been cancelled. 
Claims 241 and 246 have been objected to. 

The rejection of claims 234-247 and 466-479 have been maintained. 

The Applicants traverse the rejections and request consideration. 

As noted in the response filed March 18, 2003, all the pending rejections are based on 
ILLIADAS, either standing alone or in combination with other references. As evidenced by the 
database report from PubMed submitted along with the response filed March 18, 2003, the 
ILLIADAS reference, though having a "received" date of January 2000, has a publication date 
of June 2000. 

The Applicants respectfully attach a Declaration under 37 C.F.R. § 1.131, executed by all 
the inventors, swearing behind the ILLIADAS reference. In order to avoid duplication, the 
Applicants refer to the exhibits filed along with the unentered Declaration of Dr. Zvia Agur filed 
March 18, 2003. Specifically items 11 and 14 have been referred to. The Declaration and the 
Exhibits referred to demonstrate a conception date of at least prior to June 2000 and continued 
diligence until October 19, 2000, the filing date of the present application. 

As noted in the Declaration, the applicants completed the invention prior to June 2000, 
the date of publication of the ILLIADAS reference. Prior to June 2000, the Applicants conceived 
the invention and conveyed the subject matter to the drafters of the relevant documents attached 
to Dr. Agur's Declaration of March 2003. 
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Referring to the documents that were attached to the March 2003 Declaration of Dr. Zvia 
Agur, the Examiner incorrectly notes that the Applicants have not proven an earlier conception of 
the invention, at least that the Applicants have not proven the conception of several components 
of the invention. 

The Applicants respectfully submit that the Examiner has not reviewed all the documents 
attached to the March 2003 declaration. The Applicants attach a claim chart of the pending 
claims specifically pointing out where in the documents submitted the requisite proof of 
conception can be found. It should be noted that this is only exemplary in nature and additional 
proofs can be found in several parts of the submitted documents in addition to the ones listed in 
the attached claim chart. For convenience, additional copies of items 1 1 and 14 have been 
attached herein with page numbers for a better understanding of the attached claim chart. It 
should also be noted that the designation x:y indicates page no. x and line y in the document 
referred to. 

Therefore, the Applicants respectfully submit that ILLIADAS does not qualify as prior 
art under any sections of 35 U.S.C. § 102. 

Considering the above, it is respectfully requested that the pending claims be passed into 
allowance and the Applicants be awarded what is believed to be their rightful due. 
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CONCLUSION 

In view of the above, reconsideration and allowance of this application are now believed 
to be in order, and such actions are hereby solicited. If any points remain in issue which the 
Examiner feels may be best resolved through a personal or telephone interview, the Examiner is 
kindly requested to contact the undersigned at the telephone number listed below. 

The USPTO is directed and authorized to charge all required fees, except for the Issue 
Fee and the Publication Fee, to Deposit Account No. 19-4880. Please also credit any 
overpayments to said Deposit Account. 



Respectfully submitted, 



SUGHRUE MION, PLLC 
Telephone: (202) 293-7060 
Facsimile: (202) 293-7860 



Chid S. Iyer 
Registration No. 43,355 



WASHINGTON OFFICE 



23373 



CUSTOMER NUMBER 



Date: September 3, 2003 
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Claims 234-247 and 466-479 

A computer system for 
recommending an optimal treatment 
protocol for treating cancer using 
drugs, for an individual, said system 
interfacing with the computer and said 
system further comprising: 

a cancer system model; 



II 



III 



IV 



235. 



I 



II 



236. 



a treatment protocol generator for 
generating a_plurality of treatment 
protocols for treating cancer using 
chemotherapy; 

a system model modifier, wherein said 
the system model modifier is adapted 
to modify said cancer system model 
based on parameters specific to the 
individual; and 

a selector adapted to select an optimal 
treatment protocol from said plurality 
of treatment protocols based on the 
modified system model. 

The system of claim 234 wherein the 
system model further comprises: 



Item 11, page 70, lines 1-5 
(hereinafter 70:1-5) 



Item 11: 5:20-25, 58:20-59:10, 
54:4-7 

Item 14:2:33-40 

Item 11: 6:6-10, 64: 19-65:15, 
72:15 



Item 11: 6:1-3, 59:23-60:12, 
63:7-14,71:10-13, 72:5-7 



Item 11:6:15-17, 66:7-68:25, 
71:14-15 



a process model of cancer 
development; and 

a treatment model that is adapted to 
model the effects of treating cancer 
with drugs, including chemotherapy. 

The system of claim 235 wherein said 
process model incorporates a 
distribution of cycling cells and 
quiescent cells. 



Item 14: 13-15 



Item 11:71:7-9 



Item 14: 3:14-16 



1 



t 
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1 237. 


The system of claim 235 where a 1 Item 14: 3:18-35 
tumor cell cycle is divided into at least 
four compartments Gl, S, G2 and M 
and a quiescent stage is denoted by 
GO, wherein each of said four 
compartments is further subdivided 
into sub-compartments and an ith sub- 
compartment representing cells of age 
I in the corresponding compartment, 
wherein the system is adapted to 
ensure that cells entering a 
compartment always enter a first sub- 
compartment of the compartment. 


238. 


The system of claim 237 wherein the Item 14: 4: 15-30 
model development of cancer cells 
using a predetermined set of 
parameters by calculating a number of 
cells in each subcompartment using 
stepwise equations. 


239. 
240. 


The system of claim 238 wherein the Item 14: 3:30-35, 5:1-3 
system is adapted to use a probability 
vector is used to determine a fraction 
of cells that leaves any 
subcompartment in a compartment to 

move to a first subcompartment of the 
next compartment. 

The system of claim 238 where the Item 14- 5*4-1? ~ 

system includes a set control functions 

that are adapted to uniquely determine 

an outcome of every single step, 

wherein said control functions 

comprise age of cells, state of a 

current population and associated 

environment. 


241. 


The system of claim 238 wherein the Item 14:5:1 8-20 
system comprises a model 
representing a tumor the model 
comprising a combination of a | 
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plurality of homogeneous groups of 
cells, each of said homogeneous 
groups of cells representing a similarly 
behaving group of cells distributed 
between all the cell-cycle 
compartments. 



242, 



The system of claim 241, wherein the 
system is adapted to calculatejn each 
step, a number of cells in each sub- 
compartment of each compartment of 
each group according to factors 
including a previous state, parameters 
of tumor, tumor current 
microenvironment and drug 
concentration. 



Item 14:5:27-32 



243. 



The system of claim 242 where spatial 
structure of the tumor is included in 
the model. 



Item 14:5:32 



244. 



The system of claim 243, wherein the 
system is adapted to incorporate 
pharmacokinetic and 
pharmacodynamic, cytostatic effects, 
cytotoxic effects, and other effects on 
cell disintegration of anticancer drugs. 



Item 14:6:4-20 



245. 



The system of claim 244 wherein the 
system is adapted to incorporate^ 
dose-limiting toxicity into the model. 



Item 14:6:4-20, 7:4-20 



246. 



The system of claim 234 wherein, 
said parameters specific to the 
individual comprise parameters related 
to tumor dynamics, patient specific 
drug, pharmacodynamic and dynamics 
of dose-limiting host tissues. 



Item 14:6:4-20, 7:4-20 



247. 



The system of claim 246, wherein said 
parameters related to tumor dynamics 
comprise age, weight, gender, 
percentage of limiting healthy cells, 



Item 14:6:4-20, 7:4-20 
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desired length of treatment protocol, 
previous reaction to treatment, 
molecular markers, genetic markers, 
pathologic specifics and cytologic 
specifics. 



466. 



A computer-implemented method for 
recommending an optimal treatment 
protocol for treating cancer using 
drugs , including chemotherapy, for an 
individual, said method comprising: 



Item 11, page 70, lines 1-5 
(hereinafter 70:1-5) 



I 



creating a cancer system model; 



Item 11: 5:20-25, 58:20-59:10, 
54:4-7 

Item 14:2:33-40 



II 



enumerating a plurality of treatment 
protocols for treating cancer using 
drugs; 



Item 11: 6:6-10, 64: 19-65:15, 
72:15 



III 



modifying the system model based on 
parameters specific to the individual; 



Item 11: 6:1-3, 59:23-60:12, 
63:7-14,71:10-13, 72:5-7 



IV 



V 



467. 



selecting an optimal treatment 
protocol from said plurality of 
treatment protocols based on the 
modified system model; and 



Item 11:6:15-17, 66:7-68:25, 
71:14-15 



recommending said optimal 
treatment. 



Item 11:6:15-17, 66:7-68:25, 
71:14-15 



The method of claim 466 wherein the 
system model further comprises: 



I 



a process model of cancer 
development; and 



Item 14: 13-15 



II 



a treatment model that models the 
effects of treating cancer with drugs, 
including chemotherapy. 



Item 11:71:7-9 



468. 



The method of claim 467 wherein said 
process model incorporates a 
distribution of cycling cells and 



Item 14: 3:14-16 
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469. 



quiescent cells. 



The method of claim 467 where a 
tumor cell cycle is divided into at least 
four compartments Gl, S, G2 and M 
and a quiescent stage is denoted by 
GO, wherein each of said four 
compartments is further subdivided 
into sub-compartments and an ith sub- 
compartment representing cells of age 
I in the corresponding compartment, 
wherein cells entering a compartment 
always enter a first sub-compartment 
of the compartment. 



Item 14: 3:18-35 



470. 



The method of claim 469 wherein the 
model traces development of cancer 
cells using a predetermined set of 
parameters by calculating a number of 
cells in each subcompartment using 
stepwise equations. 



Item 14:4:15-30 



471. 



The method of claim 470 wherein a 
probability vector is used to determine 
a fraction of cells that leaves any 
subcompartment in a compartment to 
move to a first subcompartment of the 
next compartment. 



Item 14:3:30-35,5:1-3 



472. 



The method of claim 470 where a set 
control functions uniquely determine 
an outcome of every single step, 
wherein said control functions 
comprise age of cells, state of a 
current population and associated 
environment. 



Item 14: 5:4-12 



473. 



The method of claim 470 wherein a 
tumor is modelled as a combination of 
a plurality of homogeneous groups of 
cells, each of said homogeneous 
groups of cells representing a similarly 
behaving group of cells distributed 



Item 14: 5:18-20 
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between all the cell-cycle 
compartments. 



474. 



The method of claim 473, wherein in 
each step, a number of cells in each 
sub-compartment of each 
compartment of each group is 
calculated according to factors 
including a previous state, parameters 
of tumor, tumor current 
microenvironment and drug 
concentration. 



Item 14:5:27-32 



475. 



The method of claim 474 where 
spatial structure of the tumor is 
included in the model. 



Item 14:5:32 



476. 



The method of claim 475, wherein 
pharmacokinetic and [PD], cytotoxic 
effects, cytostatic effects and other 
effects on cell disintegration of 
anticancer drugs are incorporated into 
the model. 



Item 14:6:4-20 



477. 



The method of claim 476 wherein a 
dose-limiting toxicity is incorporated 
into the model. 



Item 14:6:4-20, 7:4-20 



478. 



The method of claim 466 wherein, 
said parameters specific to the 
individual comprise parameters related 
to tumor dynamics, patient specific 
drug, pharmacodynamic, and 
dynamics of dose-limiting host tissues. 



Item 14:6:4-20, 7:4-20 



479. 



The method of claim 478, wherein 
said parameters related to tumor 
dynamics comprise age, weight, 
gender, percentage of limiting healthy 
cells, desired length of treatment 
protocol, previous reaction to 
treatment, molecular markers, genetic 
markers, pathologic specifics and 



Item 14:6:4-20, 7:4-20 
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cytologic specifics. 



PATENT APPLICATION 

IN THE UNITED STATES PATENT AND TRADEMARK OFFICE 

In re application of Docket Nq . Q6068g 

Zvia AGUR, et al. 

Appln. No, 09/691,053 Group ^ Unit: 1631 

Confirmation No, 5359 Examiner: Maijorie A MORAN 

Filed: October 19, 2000 

For: SYSTEM AND METHOD FOR OPTIMIZED DRUG DELIVERY AND 
PROGRESSION OF DISEASED AND NORMAL CELLS 

DECLARATION UNDER 37 C.F.R. S 1.131 

Honorable Commissioner of 
Patents and Trademarks 
Washington, D.C. 20231 

Sir: 

We, Zvia Agur, Sarel Fleishmann, Kirill Skomorovski and Moshe Vardi hereby 
declare as follows: 

We are inventors and applicants of the invention entitled "System and Method for 
Optimized Drug Delivery and Progression of Diseased and Normal Cells", disclosed and 
claimed in U.S. Application No. 09/691,053 filed October 19, 2000. 

Prior to June 2000, we had invented the invention as described and claimed in the 
above-identified application, and pursued the present invention in the ordinary course of 
business. 

The entire inventive work was done in Israel. 

Prior to June 2000, we met with and conveyed the entire subject matter of the 
invention to the drafters of the documents included in Appendix A of the 1.131 declaration by 
Dr. Zvia Agur submitted to the USPTO on March 18, 2003. Specifically, the information 
conveyed by us prior to June 2000, also includes those appearing in items 1 1 and 14 of 



Appendix A of the 1.131 declaration by Dr. Zvia Agur submitted to USPTO on March 1 
2003. 
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The above-identified application was prepared based on the documents included in 
Appendix A of the 1.131 declaration by Dr. Zvia Agur (submitted on March 18, 2003). The 
above-identified application was then filed on October 19, 2000. 

In view of the discussion above, we had invented the claimed invention disclosed in 
the subject application prior to June 2000, and pursued the present invention in the ordinary 
course of business, including preparing the above-identified patent application, until the 
application was filed on October 1 9, 2000. 

We declare further that all statements made herein are of our own knowledge and are 
true and that all statements made on information and belief are believed to be true; and 
further that these statements were made with the knowledge that willful false statements and 
the like so made are punishable by fine or imprisonment, or both, under Section 1001 of Title 
18 of the United States Code, and that such willful false statements may jeopardize the 
validity of the application or any patent issuing thereon. 
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FIELD OF THE INVENTION 

The present invention relates to a system and method for predicting 
and optimization of treatment of disease and of certain stages in drug 
development and trials. 

BACKGROUND OF THE INVENTION 

Drugs that are administered to combat disease are often toxic to 
healthy cells as well. In prescribing a treatment protocol, it is necessary to 
consider the effect of the protocol on both abnormal and normal cells. 

Mathematical models of biological systems are well known in the art. 
As a broad category, these models are... 

More specifically, can have models of cell lines, tumor growth, etc. 
Use of these models for prediction of treatment results is also known in the art. 
However, predictive models generally employ an analytical approach, in which 
generalizations about the use of the treatment on a disease must be made. 
This approach, while providing useful general information, cannot be used to 
predict results of treatment in realistic circumstances. Thus, a method which is 
inclusive of more complex and detailed scenarios is needed. 



SUMMARY OF THE INVENTION 

To be completed when claims are finalized. 



BRIEF DESCRIPTION OF THE DRAWINGS 

The present invention will be understood and appreciated more fully 
from the following detailed description taken in conjunction with the appended 
drawings in which: 

Fig. 1 is a schematic illustration of the basis of the present invention; 
Fig. 2 is a flow chart illustration of steps of the invention, useful in 
understanding Fig. 1 ; 

Fig. 3 is a schematic illustration of a biological model, in accordance 
with one embodiment of the present invention; 

Fig. 4 is a chart illustration of the biological model of Fig. 3; 

Fig. 5 is a graphical illustration of the chart of Fig. 4; 

Fig. 6 is a chart illustration of the biological model of Fig. 3 in a 
different format; 

Fig. 7 is a graphical illustration of the chart of Fig. 6; 

Figs. 8A and 8B are graphical illustrations of the output of the model 

of Fig. 3; 

Figs. 9A and 9B are graphical illustrations of experimental data as 
compared to the output shown in Figs. 8A and 8B; 

Fig. 10 is a schematic illustration of a biological model, in accordance 
with a further embodiment of the present invention; 

Fig. 1 1 is a graphical illustration of results of the simulation of the 
model shown in Fig. 1 0; 

Figs. 12A and 12B are graphical illustrations of the effects of two 
doses of G-CSF on the neutrophil line, according to the model of Fig. 10; 



Fig. 13 is a schematic illustration of a biological model, in accordance with 
further embodiment of the present invention; 



DETAILED DESCRIPTION OF THE PRESENT INVENTION 

A system and method have been developed for identifying optimal 
treatment protocols for selected parameters. Unlike prior art optimization 
schemes, the present invention uses heuristic determinations. The use of 
heuristics enables to find near optimal solutions even for complex (hence 
realistic) mathematical descriptions of the relevant 
biological/clinical/pharmacological scenarios; this holds both for the general 
case ("optimal generic treatmenf), as well as at the level of an individual 
patient. Reference is now made to Fig. 2, which is a flow chart illustration 
generally showing the optimization method. The Figure generally depicts the 
basic concept. The present invention optimizes a drug delivery protocol after 
consideration of possible protocols, referred to as a protocol space. 
Determination of optimal protocol is based on specific parameters input by a 
user . The user may be a physician, a drug developer, a scientist, or anyone 
else who may need to determine a treatment protocol for a drug. The specific 
parameters may include several categories, such as individual patient 
characteristics and/or medical history, needs of a specific user (research, 
efficacy, treatment, etc.), and other particulars (such as maximum length of 
treatment, confidence level, etc.). 

The first steps include generation of several types of models. A set of 
models for all the relevant biological processes, is created. In addition, a 
model of treatment effects on each of these processes is created. The 
combination of these models provides a detailed mathematical model of the 
overall bio-clinical scenario in a general or specific sense together with the 
specific effects of a particular treatment. Once the comprehensive model is 



constructed the characteristic parameters (either population averaged or 
patient specific e.g., - age, gender, weight or clinical indications) are 
implemented in it. In this way a 'Virtual" patient is generated.The next step is 
definition of a protocol space. To do this, possible permutations of certain 
parameters such as drug doses, dosing intervals, drug quantity, etc. are 
considered. Thus, a number (can be very large) of possible treatment 
protocols is generated. The amount of possibilities depends on the number 
and ranges of parameters considered. At this point, the fitness function is 
constructed by mathematically considering different possible factors which 
may be influenced by the treatment. These may include survival, tumor or 
pathogen load, cytotoxicity, side effects, pain etc. The user can put certain 
specific parametersin the fitness function so as to adjust this function to 
his/her specific goalsBased on the selected parameters, the fitness function is 
applied, and calculates a fitness score for each and every possibility in the 
protocol space. Finally, the optimization step is carried out, either by search 
heuristics or by analytical methods, in order to select the optimal treatment 
protocol from all the scored possibilities. 

In this way, a disease specific, patient specific, situation specific, drug 
specific, objective specific treatment protocol may be obtained. The actual 
time it takes once the parameters are entered may be negligibly short or up to 
hours, depending on the length of the simulated treatment period and the 
power of the specific search heuristics and the computational tools, making 
this a very feasible tool. 

Construction of detailed mathematical models for biological processes 
and treatments will be shown by way of example. A model of platlets 



production and thrombocytopeniawith treatment by TPO, a model of neutrophil 
production and neutropenia and treatment by G-CSF, and a model of cancer 
growth and chemotherapy are described. An example for specific optimization 
(by linear programming) is implemented in the neutrophil model. A general 
heuristic optimization method is described as well. 

EXAMPLE ^Thrombocytopenia and Treatment by Thrombopoietin (TPO) 

Thrombocytopenia is a common hazardous blood condition, which may 
appear in different clinical situations, including cancer chemotherapy. 
Recently, a thrombopoiesis-controlling cytokine, thrombopoietin (TPO), was 
isolated and its human recombinant analog became available. A 
mathematical model has been developed to imitate dynamics of a 
thrombopoietic lineage in the bone marrow, of platelet counts in the periphery, 
and effects of TPO administration on them [on what? Platelets?]. 

TPO is a cytokine, glycoprotein of about 350 amino acids, that 
resembles erythropoiesis-stimulating hormone, erythropoietin. Its synthetic 
analogs, recombinant human thrombopoietin (rHuTPO) and recombinant 
human megakaryocyte growth and development factor (rHuMGDF), are 
available as well and are undergoing clinical trials. These compounds have 
been shown to have the same biological activity as TPO has, so the term 
TPO will be used without distinguishing between its different forms and 
analogs. 

TPO is a primary growth factor of the thrombopoietic cell line both in vivo 
and in vitro. Aside from this, TPO may be potent in stimulation and co- 



stimulation of other hemopoietic lineages (e.g., granulopoietic or 
erythropoietic). 

I. MODEL OF BIOLOGICAL SYSTEM 
A Background of Thrombopoiesis 

Like all other hemopoietic lines, the thrombopoietic line originates from 
poorly differentiated, multipotential cells, capable of some division and self- 
reconstitution. Other bone marrow cell lines, such as pluripotential 
hemopoietic stem cells and common myeloid progenitor cells (CFU-GEMM) 
have these characteristics as well. 

Starting out as undifferentiated stem cells, the cells gradually become 
more and more differentiated and committed to a specific line. Once they 
commit to the thrombopoiesis line, they proliferate extensively within certain 
compartments. One such compartment is that of colony-forming units - 
megakaryocytes (CFU-Meg). Other examples of compartments burst-forming 
units - megakaryocyte (BFU-Meg) and megakaryoblast compartments, which 
have similar properties as one another. 

The committed thrombopoietic cells, which are megakaryocyte 
precursors, go through several stages of maturation when proliferation is no 
longer possible. However, megakaryocyte maturation is somewhat different 
from maturation within other hemopoietic lines. In megakaryocytes, cell nuclei 
undergo mitosis in parallel with cytoplasmic maturation. However, although 
the DNA material of these cells doubles, cell division does not happen. Such 
incomplete mitosis is called endomitosis. Consequently, the cell becomes 



0 

polyploid with 2N, 4N, 8N, etc., amounts of DNA. The cells with 2N to 4N 
chromosomes are called either megakaryoblasts or megakaryocytes. 

Usually, megakaryocytes do not start to release platelets until the 16N 
state. Then they begin to create demarcation membranes that envelop 
5 cytoplasm fragments generating platelets [what does this refer to?]. The 
platelets are then released into the blood stream. The amount of cytoplasm, 
cell volume and an ability to release platelets all increase proportionally to the 
cell ploidy. 

B. Mathematical Model 

10 Reference is now made to Fig. 3, which is a detailed illustration of a 

model predicting thrombopoiesis. As shown in Fig. 3, the thrombopoietic 
lineage is divided into eight compartments. The first compartment, called 
Stem Cells (SC) and labeled 30, refers to all bone marrow hemopoietic 
progenitors that have an ability to differentiate into more than one line (e.g., 

15 pluripotential stem cells, CFU-GEMM, and others). Cells of SC compartment 
30 proliferate, mature, and subsequently differentiate into megakaryocytes or 
other precursors, or they may give rise to new stem cells. Although the term 
"new" is not biologically accurate, it serves as an acceptable assumption for 
the purposes of the model. 

20 Cell death through apoptosis may have a significant effect on cell 

number within proliferating compartments. However, the effect of apoptosis is 
combined with the effect of cell proliferation into a total amplification of cell 
number in a given compartment. Thus, for example, the amplification rate of 
cells in SC compartment 30 is normally 1.029 new cells per hour. An 



1 
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assumption is made that no apoptosis occurs in non-proliferating 
megakaryocyte compartments, due to lack of evidence to the contrary. 

Biologically, rates of proliferation and maturation, the ability to 
reconstitute, and other characteristics differ between particular cell types 
5 within a primitive progenitor population. However, in this model there is no 
distinction between them; all progenitor cells are considered to be one 
population with common properties. 

Based on some works [which ones?] that have shown that probabilities 
of stem cell differentiation into one or another hemopoietic lineage are 
10 constant, it is assumed that a flow of stem cells into the megakaryocyte 
lineage depends only on the number of mature stem cells. The same was 
assumed about the stem cell self-renewal. Thus, after the cells spend a 
defined transit time, for example 24 hours, in SC compartmen^30) a certain 
constant fraction of the cells (0.5181) return to their "young state", i.e. start 



15 their passage through SC compartment 30 again, /aSshown in line 31 



Another constant fraction (0.342, for example) of cells pass into the next 
compartment called Colony-Forming Units (CFU-Meg), labeled 40. It is 
presumed that remaining stem cells differentiate into other hematopoietic 
lineages. 



megakaryocyte line but are still capable of proliferation. Cells of CFU-Meg 
compartment 40, like those of SC compartment 30, spend some time 
(normally 60 hours) multiplying at an amplification rate of about 1.0255 [per 
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CFU-Meg refers to all cells that are already committed to the 



It) 



what?] and maturing before losing their proliferative abilities and passing on 
to the next compartment 50, called megakaryoblasts (MKB). 

MKB compartment 50 includes all the cells that have lost the ability to 
proliferate, but are not yet sufficiently mature to release platelets. For the 

5 purposes of the model, the assumption is made that megakaryocytes do not 
start to release platelets until they reach the 16N-ploidy phase. Hence, MKB 
refers to 2N, 4N and 8N cells of megakaryocyte lineage that cannot divide, at 
all stages of cytoplasmic maturity. After these cells spend the designated 
transit time (normally 132 hours) in MKB compartment 50, they move to the 

10 next compartment 60, which is an MK16 bone marrow compartment. 

The cells of MK16 compartment 60 are megakaryocytes of 16N-ploidy 
class that release platelets until they exhaust their capacity, and then are 
disintegrated. Cell volume has a linear relationship with megakaryocyte 
ploidy. Hence, we assume that all 1 6N-megakaryocytes have the same 
is volume and, thus, the same platelet-releasing capacity. Therefore all platelet- 
releasing 16N-megakaryocytes are in transit for the same amount of time 
(120 hours) until they are exhausted and disintegrated. 

However, some 1 6N-megakaryocytes do not participate in platelet 
release, but rather continue with another endomitosis over a 48-hour time 
20 period, and become 32N-megakaryocytes. These constitute a new and 
distinct MK32 compartment 70. Thus, after 48 hours in MK16 compartment 
60, a certain fraction of the cells (normally 0.26) leave MK16 compartment 60 
and go on to MK32 compartment 70. 



32N-megakaryocytes release platelets as well. The rate of platelet 
release is constant for every compartment and proportional to the ploidy state 
of megakaryocytes in it. Thus, every 1 6N-megakaryocyte releases 10.2 
platelets per hour and every 32N-megakaryocyte releases 20.4 platelets per 
hour. However, 32N-megakaryocytes are not exhausted more quickly than 
16N-megakaryocytes, since they have 2 times greater volume and platelet- 
releasing capacity. Consequently, all platelet-releasing megakaryocyte 
compartments have the same transit time. 

Once again, some fraction of cells (normally 0.016) are not engaged in 
platelet formation, and continue to the 64N-stage. Additional endomitosis in 
MK32 compartment 70 takes the same amount of time (48 hours) as in MK16 
compartment 60. The 64N-megakaryocytes continue the process in a new 
MK64 compartment 80, and 0.4 of them become 128N-cells in yet another 
MK128 compartment 90. Additional endomitosis in MK64 compartment 80 
takes the same amount of time (48 hours) as before. 

Normally, 64N- and 128N-megakaryocytes are not found in humans. 
However, these cells were detected in human bone marrow in some abnormal 
conditions. Therefore, we included these compartments in our model, but in 
the normal state their number is negligible. As far as we know, 
megakaryocytes of greater ploidy classes have not been encountered in 
humans. 

Finally, there is a platelet compartment 100. This is not a bone marrow 
compartment, but rather the platelet pool in the peripheral blood [(spleen 
sequestration???)]. Platelets released from megakaryocytes of 16N-, 32N-, 



64N-, and 128N-ploidy classes enter platelet compartment 100. Platelet 
elimination from the blood is presumed to occur after a defined transit time in 
circulation by a mechanism of death from senescent [what is this?]. To 
simplify the model we ignore the destruction of platelets by their use [?]. 
Transit time through platelet compartment 100 is 240 hours but can vary, as 
will be discussed below. 

II. MODEL OF TREATMENT EFFECTS 
A. Background of TPO 

The major sites of TPO production are the liver and kidney. TPO is also 
produced in the spleen and bone marrow, but the production rate in these 
organs is 5 times lower than in the liver and kidney. Some low TPO 
production has also been found in many other sites in the body. Rates of liver 
and kidney TPO production are constant under thrombocytopenia and 
thrombocytosis of varying degrees of severity. TPO production in the spleen 
and bone marrow is inversely related to the megakaryocyte mass, but the 
actual contribution is negligible with regard to total TPO production. 

Another mechanism of TPO concentration regulation is receptor- 
mediated TPO uptake, since TPO-receptors on the platelet and 
megakaryocyte surfaces are the main TPO-clearance mechanism. Thus, 
TPO concentration is inversely related to the total platelet and megakaryocyte 
mass. 

The effects of TPO on the thrombopoietic line may be divided into three 
types: (i) stimulation of proliferation of megakaryocyte progenitors that have 



an ability to proliferate; (ii) stimulation of maturation of all megakaryocyte 
progenitors; (iii) induction of additional endomitosis of already mature 
megakaryocytes, which leads to an increase in the modal megakaryocyte 
ploidy. 

B. Mathematical Model of TPO Effects 

TPO concentration effects on the thrombopoiesis line is now considered. 
As discussed above, three things depend on TPO concentration: (i) 
amplification rate (amp), (ii) the rate of cell maturation or, alternatively, transit 
time through a given compartment (transit), and (iii) the fraction of 
megakaryocytes of given ploidy that undergo additional endomitosis and pass 
on to the next ploidy class. 

1 . TPO Concentration 

Recombinant human full-length TPO and its truncated form rHuMGDF 
are fully active biologically. Therefore, in our model we add exogenously 
administered recombinant protein to endogenously produced TPO in order to 
calculate actual TPO concentration. 

As mentioned above, the rate of TPO production in the main TPO 
production sites, i.e. liver and kidney, is constant under thrombocytopenia or 
thrombocytosis. The level of TPO mRNA in sites like the bone marrow and 
spleen, where it is produced in a 5-fold lower rate than in the liver and kidney, 
is not significantly different from the TPO level in peripheral blood. Therefore, 
the assumption is made that the bone marrow and spleen contributions to the 
total TPO concentration are insignificant. Endogenously produced TPO is 



assumed to have a constant rate of production of 7 pg/ml/hour. However, this 
number can change. 

The main mechanism that controls TPO concentration in the blood is 
5 receptor-mediated TPO uptake. Both megakaryocyte and platelet mass 
contribute to the total receptor number (normally 8.375 x 10 9 ) and, 
consequently, to the rate of TPO clearance. We assume that every platelet 
bears 220 TPO receptors [based on what?]. Each megakaryocyte at any 
given moment bears an amount of receptors equal to the number of platelets 
10 that can be released during the remaining life-span, times 220 receptors per 
potential platelet. We assume also that every receptor molecule removes 
4.776 x 10~ 10 pg/ml/hour of TPO from circulation. However, this number can 
change. 

Another mechanism of TPO removal from the blood is non-specific 
15 TPO- receptor-independent clearance. This mechanism is rather insignificant 
in the normal state, when receptor-mediated TPO binding, endocytosis, and 
degradation remove most of the TPO. However, when the amount of TPO 
rises significantly above the ability of the receptor pool to remove it, non- 
specific clearance becomes important. In our model this type of TPO 
20 clearance is exponential, i.e. every hour some fraction (0.03) of a current 
amount of TPO is removed from circulation. 

Exogenous TPO is included in the model as a linear relation of the 
initial maximum TPO blood concentration to the administered intravenous (IV) 
dose (the relation coefficient is 0.0167) based on literature data. 



The state when TPO completely disappears from the blood is very 
unlikely, the lower limit of possible TPO concentration is restricted to a certain 
minimum, in this case 0.01 pg/ml. 

Thus, the formula that calculates TPO concentration hourly is given in 
Equation 1 as follows: 

C M = max((C;. + Endo+ 0.0 1 67- Exo- Nc • C, - Rc • /?p),0.0 1) (1 ) 

where C, is TPO concentration at the current hour (/), C i+1 is TPO 
concentration at the next hour, Endo is endogenously produced TPO per 
hour, Exo is the total amount of TPO administered intravenously during the 
current hour, Nc is the coefficient of non-specific TPO clearance, Rc is the 
coefficient of receptor-mediated TPO clearance, i.e. amount of TPO that each 
receptor removes per hour, and Rp is the total TPO-receptor pool. 

In steady-state TPO concentration is constant and equals 100 pg/ml. 

2. TPO effects on amplification rate 

In our model, there are only two compartments, SC compartment 30 and 
CFU-Meg compartment 40, whose cells are capable of dividing. These 
compartments differ significantly from each other, so they will be discussed 
separately. Cells of other compartments do not proliferate, so their 
amplification rate equals 1 under all circumstances. 

a. SC compartment 30: 

Cells of SC compartment 30 are relatively insensitive to TPO as 
compared to committed megakaryocyte cells. In our model this is considered 



as a threshold of 2500 pg/ml of TPO concentration, but this number can be 
changed. TPO affects stem cells only above this concentration. As long as 
TPO remains below the threshold, stem cells are regulated by an intrinsic 
TPO-independent mechanism that keeps the size of their population almost 
5 constant. 

Thus, below the threshold, SC amplification rate (ampsc) is determined 
hourly depending on the current number of cells in SC compartment 30. The 
dependence equation gives a sigmoidal function: ampsc changes from 1 (i.e., 
no amplification, the cell number remains the same) up to the maximal value 
10 when the cell number approaches infinity or 0, respectively [??]. When the 
cell number is normal [define normal] amp S c equals one fourth of its 
maximal possible value. Stem cell amplification rate is defined by Equation 2 
as follows: 



15 / . \ ( ^ quantnorm sc s (2 ) 

amp sc [quant sc ) = [amp SCmM - 1)-- 5 j- + 1 

3 quant sc + quantnorm sc 



where ampscmax is the maximal possible rate of cell amplification in SC 
compartment 30, quantnorm S c is the normal quantity of cells in SC 

20 compartment 30, quantsc is an actual quantity of cells there, and S is a value 
that determines the sensitivity of the mechanism that links the amplification 
rate to the number of cells. A high value for S means that ampsc changes 
significantly due to silent changes of quantsc, and a low value for S means 
that ampsc remains relatively constant whatever the changes of quantsc may 

25 be. S generally equals 0.2, but this number can be changed. 
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In a range of very low cell numbers (100 times lower then normal or 
less), the model must define the stem cell parameters for the quickest 
recovery of SC compartment 30. Thus, amp S c in these conditions achieves 
the maximal value (ampscmax). 

5 On the other hand, the rise of TPO concentration above the threshold 

should occur in severe platelet and/or megakaryocyte deficiency, or when 
TPO has been administered exogenously. In these situations the hemopoietic 
system mobilizes additional resources for recovery of platelets and platelet- 
releasing compartments. Hence, we assume that in an effective 

10 concentration, TPO increases the rate of cell amplification in SC compartment 
30. This is shown in Equation 3, as follows: 



amPsc(<l uant sc* C i) = 



/r / v quantnorm sc s 



3 • quant sc s + quantnorm sc s 



+ 1 



(amp SCn ^ -l).F*(G # (C, ~thres\St) 



(3) 
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where the first of two operands is the amplification calculated based on a 
TPO-independent mechanism, and the second is the TPO- related addition to 
the amplification rate. 

G* is a function that transforms concentration (C) into a form that is 
20 easy to use in the calculation of both ampsc and transitsc as will be discussed 
futher [where?]. G* is defined by Equation 4 as follows: 



G* (c) = ln(C + 1) • 0.23 + 0.03949 



(4) 



25 



F* in Equation 3 is a function that determines the steepness, denoted Sf, of 
the ampsc-versus-G curve. F, as shown in Equation 5, is a recurrent function 



that adds 1 to its first argument and applies a logarithm to this sum. 
Thereafter, 1 is added another time and a logarithm is applied to the sum 
again. This operation is repeated Sf times. 

F * (x ) = ... log (log (log (log (* + 1 )+ 1 )+ 1 )+ 1 ).. (5) 

Although St appears in several equations, its value is specific for every case. 
In Equation 3, thres is the threshold above which TPO has an effect on the 
cells. In order to avoid an abrupt change of amp sc when TPO concentration 
transcends the threshold, the TPO-related addition is calculated relative to the 
difference between an actual TPO concentration and the threshold. 

b. CFU-Meg Compartment 40: 

In contrast to the cells of SC compartment 30, cells of CFU-Meg 
compartment 40 are very sensitive to TPO and respond to the absolute TPO 
concentration, rather than its level above a certain threshold. Biologically, 
these cells have no TPO-independent proliferation mechanism, and thus 
cease to proliferate when deprived of TPO. However, in our model, CFU-Meg 
cells continue to proliferate without TPO at almost the same rate as with TPO. 
With normal TPO concentrations, ampcFu equals one eighth of its maximal 
possible value (ampcFUma*)- The equation that describes the relation of 
amplification rate of CFU-Meg cells (amp C Fu) to TPO concentration (CI) is: 

amp CFU (C, )= (amp CFU max - l)- F(G(C, \St)+l m 



Function G resembles function G* (Equation 4) qualitatively but differs 
quantitatively, having the following form: 

G(c)=2.56xl0" 6 -(ln(C + l)+ 0.3949 ) 7 +0.3 (7) 

Function G is changed from function G* by adding 0.3949, raising to a power 
of 7, multiplying by 2.56x1 0' 6 , and adding 0.3. These changes provide the 
following features: 

(a) in the region near the normal TPO concentrations the function 
rises with an increasing rate from about 0.3 up to about 0.8; 

(b) in the normal TPO concentration (100 pg/ml) it equals 0.5; 

(c) in high TPO concentrations, which occur in cases of exogenous 
TPO administration, the function rises with a decreasing rate; 

(d) despite (iii), in high TPO concentrations that may arise in a 
logical situation (i.e. severe thrombocytopenia or TPO 
administration in logical doses), the rate of function growth does 
not drop down to zero and the function does not reach a 
maximum, but rather the rate gradually stabilizes at a nearly 
constant value and the function continues to rise almost linearly. 

Function F in Equation 6 is almost the same function as function F* in 
Equation 5. The difference is that while F* returns the number solely on the 
basis of Equation 3, F is required to return to 0.5 when its first argument is 
0.5, no matter what St is. This is achieved by multiplying the returned number 
by a certain coefficient, which is specific for every value of St. Thus, F is 
defined by Equation 8 as follows: 



F(x)=k- (... log (log (log (log (x + l)+ l)+ l)+ 1)..) 
where k is the Sf-specif ic coefficient. 



(8) 



3. TPO effects on transit time 

As noted earlier, all platelet-releasing megakaryocyte compartments 
have the same transit time. Since neither the relationship of megakaryocyte 
volume (and thus, its platelet releasing capacity) to megakaryocyte ploidy, nor 
its rate of platelet release, are affected by TPO, the transit time through the 
compartments are constant, for example 120 hours. Platelets in circulation 
also spend constant period of time (eg. 120 hours) which is not affected by 
TPO concentration. 

In contrast, SC, CFU-Meg, and MKB compartments have changeable 
transit times. If the transit time were continuously strictly related to the TPO 
concentration, then during an abrupt concentration change (in exogenous 
administration, for example) there could be either massive cell exit from the 
given compartment or absolute absence of cell exit for the period of time that 
equals the difference between current and previous transit time. Since this 
situation seems biologically unlikely, in the present model transit time 
changes gradually, approaching the value that it should be equal to based on 
TPO concentration step by step. 

Regarding the amplification rate, transit time in SC compartment 30 
differs from the other compartments. 

a. SC compartment 30: 



Cells of this compartment respond to TPO only when its concentration 
rises above the threshold. Below the threshold stem cell transit time is 
regulated by an intrinsic mechanism dependent on the current cell number 
only. The value that the transit time should approach is determined based on 
the cell number, according to the following equation [=?]: 




(9) 



where transitscmin is the minimal possible transit time through SC 
compartment 30 (1 2 hours). It is two-fold lower than the normal transit time. 

Unless the value that transit time should approach does not differ 
from the current transit time by more than 1.0, transit time does not 
changes. When it does differ, the transit time begins to change every 
hour by 1 hour in the direction of the value, [explain??] 

When TPO concentration rises above the threshold, it sets the value 
which stem cell transit time approaches, to the following number [=?]: 




(10) 



where G* is the function given in Equation 4, C, is the actual TPO 
concentration and thres is the same threshold as for stem cell amplification 
rate. Thus, TPO shortens stem cell transit time, which allows for the quickest 
recovery of committed megakaryocyte compartments when those [what?] are 
diminished. [?] 



1 



1 



We assume maximal possible transit time (300 hours), above which 
transit time does not rise. 

As in amplification rate, in a range of very low cell numbers (100 times 
lower then normal or less), the model must define the stem cell parameters 
5 for the quickest recovery of SC compartment 30. Thus, transitsc in these 
conditions achieves the maximal value, which is 300. 



b. CFU-Meg and MKB compartments: 



Transit times of these CFU-Meg and MKB compartments 40 and 50 are 
10 completely dependent on TPO concentration. The value, which the transit 
time approaches, is [=?]: 



15 




(11) 



20 where transit (C Fu-MegMKB) mm is the minimal possible transit time through CFU- 
Meg compartment 40 (eg. 30 hours) or MKB compartment 50 (66 hours), 
respectively. G is given in Equation 7, and F is the same function as in 
Equation (7 or 8??), but the value of St is different in each equation. Maximal 
possible transit time is 100 hours for CFU-Meg compartment 40 and 350 

25 hours for MKB compartment 50. 

4. TPO effects on cell flow from one compartment to another 

Cell flow between compartments refers not to the rate of cell passage 
from one compartment to the next, and not to the number of cells that pass 



during a time unit, but rather to the proportion of "mature" cells that pass to 
the next compartment at any given moment, designated "floworf. "Mature" 
cells are ones that are potentially ready to pass to the next compartment but 
do not necessarily pass. 
5 As was noted earlier, it is assumed that the fraction of SC compartment 

30 that commits to the megakaryocyte lineage is constant (0.342) and TPO- 
independent. From the two following compartments, CFU-Meg compartment 
40 and MKB compartment 50, every mature cell emerges to the next 
compartment regardless of external circumstances [is this correct? Or just 

10 TPO circumstances?]. Thus, TPO does not affect the flowon in the first 
three compartments. 

In contrast, the fraction of MK16-, MK32-, and MK64- megakaryocytes 
that continue with additional endomitosis and flow to the next compartment 
may change from almost 0 up to 1 depending on TPO concentration. 

is Because there is no compartment with ploidy greater than 128N, the 
megakaryocytes of MK128 compartment 90 do not flow to any other 
compartment, so there is no flowon value for MK128 compartment 90. 

The dependence of MK16, MK32, and MK64 flowon parameters on TPO 
concentration is expressed by the following bi-phasic function: 



20 




(12) 



25 where flowonnorm is the value of flowon under normal TPO concentration 
(100 pg/ml). This value for MK16 is 0.26, for MK32 it is 0.016, and for MK64 
it is 0.4. 
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G° is a function that resembles the function G in Equation 7. The 
difference is that G° relates to the TPO concentration before a certain period 
of time rather than to the current TPO concentration. The reason for this is 
that it is assumed that a cell that enters a given megakaryocyte compartment 
first "decides" whether to undergo additional endomitosis and not participate 
in platelet release in this compartment, or whether to begin with platelet 
release and remain in this compartment until complete exhaustion. TPO- 
dependent determination of what fraction of cells will "choose" each possibility 
occurs at the start of the cells' path through the given compartment. However, 
the result of this determination can be seen when the cells that "decide" to 
leave the compartment actually leave it, i.e. after they complete one 
endomitosis. Thus, the fraction of cells that pass to the next compartment, 
flowon, is calculated based on the TPO concentration that existed before the 
period of one endomitosis (48 hours): 




(13) 

The time needed for additional endomitosis is the same in all three 
compartments, MK16, MK32, and MK64. 

HI. COMPLETE DETAILED MODEL 

The complete model was built as an imitation of what happens in real 
bone marrow. Each compartment is subdivided into small sections that 
contain the cells of a specific age with a resolution of one hour. For example, 
the fifth age-section of MKB compartment 50 contains cells within MKB 
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compartment 50 that have been within that compartment for 5 hours. Every 
hour, all the cells in the "bone marrow" pass to the next age-section in the 
same compartment. 

When the cell has spent all the transit time predetermined for it in a 
given compartment, it passes to the next compartment to the zero age- 
section. Thus, every hour the cells that leave one compartment fill the zero 
age-section of the next one. The cells that leave MK128 compartment 90 
have nowhere to go and therefore disappear. The zero age-section of SC 
compartment 30 is filled by a certain fraction (0.5181) of the cells that leave 
SC compartment 30. 

The cells that release platelets add a certain platelet number to the zero 
age-section of PL compartment 100 every hour. 

Reference is now made to Fig. 4, which is an illustration of the 
implementation of the model. The model is implemented as a chart of 8 rows 
and 360 columns. The 8 rows relate to 8 cell compartments, and the columns 
relate to the age sections, with the assumption that transit time does not 
exceed 360 hours. This chart is updated hourly according to the rules 
described above. 

Reference is now made to Fig. 5, which shows a graphical 
representation of the chart of Fig. 4. Within the compartments that 
proliferation occurs (SC and CFU-Meg), the number of proliferating cells 
increases from the first to the last age-section. In contrast, the cell number in 
the compartments that have no proliferating ability remains constant (MKB, 
MK128, PL), or decreases when cells that have undergone additional 
endomitosis leave the compartment for the next one (MK16, MK32, MK64). 



Reference is now made to Fig. 6, which is an illustration of another 
representation of the model, based on the time courses of different 
compartments. The rows in the chart represent cell compartments and the 
columns represent time of simulation course. At every time-step of the 
simulation (one hour of "patient's life"), the number of cells in all age-sections 
is summarized for each compartment and the next column in time-course 
chart (Fig. 6) is filled. Thus, every cell in the chart represents the total 
number of cells in a given compartment at a given time point. 

There is an additional row in the time-course chart that relates to the 
TPO concentration in the blood. TPO concentration is written down every 
time-step concurrently with the cell numbers. 

Reference is now made to Fig. 7, which is a graphical representation of 
the chart of Fig. 6, and is the most useful model output [why? Describe 
graph and how it is used]. 

The implementation of the described model results in a computer 
simulator that describes the changes that occur in the human thrombopoietic 
system (platelet counts, bone marrow precursor numbers, and TPO 
concentration) over several years. The resolution of the simulator output is 
one hour. 

Time units and periods that we will speak about relate to the 
simulated patient's life, not to the running time of the program. 
IV. PARAMETER-SPECIFIC ADAPTATION OF MODEL 

This model may be fit to patients with diverse blood and bone marrow 
parameters. People differ in their baseline platelet counts and numbers of 
bone marrow precursors, in the cells' transit times and amplification rates, 



rates of platelet release by megakaryocytes, fractions that each 
megakaryocyte ploidy class donate for additional endocytosis, and in the time 
needed for endomitosis. Furthermore, the rate of TPO production, receptor- 
and non-receptor-mediated TPO clearance, the threshold of TPO effect on 
5 SC compartment 30, and the sensitivity of different cell parameters to TPO 
also differ from one person to the next. 

To obtain an ideal fitness of the model to each patient, the patient- 
related parameters should be given individually for each patient. However, 
practically, it would be extremely difficult to predetermine many of these 

10 parameters for every patient. Therefore, certain average parameters have 
been calculated based on published data, and are shown in Table 1 below. 
These averaged parameters are used as a framework into which known 
individual characteristics are included. Thus, before a particular simulation is 
begun, relevant known information about the individual may be included, 

15 sometimes replacing certain parameters of the model. 



TABLE 



^^Feature 
Compartment^^ 


Baseline 
number 
(xl0 3 /kg) 


Normal 
amplification rate 
(increase per hour) 


Normal transit time 
(hrs) 


Normal cell 
flow* 
(hour - 1 ) 


Rate of platelet 
release (hour -1 x 
megakaryocyte _1 ) 


Time needed for 

additional 
endomitosis (hrs) 


SC 


478 


1.029 


24 


0.342 






CFU-Meg 


1250 


1 .0255 


60 


1 






MKB 


5105 




132 


1 






MK16 


4080 




120 


0.26 


10.2 


48 


MK32 


1250 




120 


0.016 


20.4 


48 


MK64 


15 




120 


0.4 


40.8 


48 


MK128 


7.5 




120 




81.6 




PL 


15 050 000 




240 








TPO-related 
parameters 


Rate of 
production 


Non-receptor- 
mediated clearance 


Receptor- 
Mediated clearance 


Threshold for 
the effect on SC 








7 pg/ml/hr 


3 %/hr 


4.776 x 10- 10 
pg/ml/receptor/hour 


2500 pg/ml 






Steepness of the 
TPO sensitivity 
curve of 


SC 
amplifica 
tion rate 


CFU-Meg 
amplification 
rate 


CFU-1S 
transit i 


rteg MKB 
time transit time 


Fraction of MK16 
undergoing 
additional 


Fraction of MK32 
undergoing 
additional 


Fraction of MK64 
undergoing 
additional 



1% 



different 
parameters (St) 










endomitosis 


endomitosis 


endomitosis 






30 


100 


7 


3 


3 


3 



* - fraction of mature cells of a given compartment that goes to the next compartment 

** - Sensitivity of the intrinsic TPO-independent mechanism that determines SC amplification rate (S) is 0.2. 

*** - the logarithm base is 10, not e, like in others 

5 



Usually, the known patient-related data are not parameters in the form 
defined by our model, but rather measurements obtained in the clinic (e.g., 
day and value of post-chemotherapy thrombocytopenia nadir, day and value 

10 of platelet peak after TPO administration, change in megakaryocyte modal 
ploidy following some perturbation, etc.). In these cases, the available data is 
converted into a model-compatible format. 

Sometimes, the only available patient-related data is the graphic 
representation of the patient's platelet course following some perturbation 

15 (e.g., cell-suppressive therapy or TPO administration). The data may also be 
a picture of the platelet course without any external disturbance (e.g., cyclic 
thrombocytopenia). In these cases the model parameters are changed by 
trial-and-error until a good compliance of the model graphic output and the 
patient's graphs is achieved. It should be noted, however, that even in the 

20 case of trial and error, the choices of parameter sets are not random but 
rather are also based on some analysis [for example?]. 

Specifically, the following tools are available for providing maximum 
flexibility: 

1) The user can set the baseline values and all other known patient- 
25 specific thrombopoietic parameters before starting the simulation. 

2) The user (e.g., physician) can determine how long of a time period to 
simulate, from 12 hours up to several years. 



3) The user can determine the frequency of showing the course of a 
patient counts up to the moment. The frequency can change from as 
much as every 12 hours to once during the overall period of 
simulation. 

4) The user can determine the resolution of the output graph, from the 
hourly representation of the patient's state down to any other 
resolution. [What does this mean? Can it be resolved at units of 
less than an hour? If so, how?] 

5) The user can choose to view the graphical representation of the age 
distribution through the compartments at any moment of the 
simulation. 

6) The user can imitate a cell-suppressive therapy at any moment while 
running the simulation by reducing one or several of the 
compartments by any value. 

7) The user can simulate exogenous TPO administration at any 
moment while running the simulation by controlling dose height, 
number of dosings or frequency of dosings. 

The simulation tool has been carefully tested with respect to the 
published experimental results, and has proved to be well calibrated for 
average human data. Parameters may be modified relatively quickly for 
efficient use of the system. The following model parameters are important for 
individualized adjustment of the model: 

• baseline number of: SC, CFU-Meg, MKB, MK16, MK32, platelets. 

• amplification rate of: SC, CFU-Meg. 
•transit time of: MKB, MK16, MK32, MK64. 



•fraction undergoing additional endomitosis in: MK16, MK32, MK64. 

• rate of platelet release of: MK16, MK32, MK64, MK128. 
•Time needed for additional endomitosis. 

• Rate of endogenous TPO production. 

5 • Ratio of receptor- and non-receptor-mediated TPO clearance. 

• Steepness of the sensitivity curve of: CFU-Meg amplification rate; MKB 
transit time; MK16, MK32, and MK64 fraction undergoing additional 
endomitosis. 

[Question: which of these parameters can actually be measured 
10 and input into the program? All? Are there other parameters too 

such as age, weight, medical history, etc.? I think we need to 
distinguish between types of parameters.] 

Reference is now made to Figs. 8A, 8B, 9A and 9B, which show a 
comparison between experimentally obtained data and the simulated model. 

15 Experimentally obtained in vivo platelet counts following TPO administration 
are shown in Fig. 8A [is this with chemotherapy too?], and chemotherapy 
without TPO is shown in Fig. 9A. Figs. 8B and 9B show simulations of the 
same. By using a TPO schedule designed by the described model, one can 
obtain platelet profiles that are similar to those obtained clinically (Fig. 8B) or 

20 even more effective (Fig. 9B). In this case, these results are achieved by 
administering a pre-calculated TPO protocol whose total dose amounts to 
25% of the original total dose, explain the point of this better 

The complete model simulates cell and platelet counts in the steady 
state, as well as after perturbations to the hematopoietic system, e.g., cell- 

25 suppressive therapy, recombinant thrombopoietin administration, etc. It is 
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possible to simulate any protocol of drug administration and any 
hematological state of a patient, regarding his/her platelet count and number 
of bone marrow megakaryocytes and their precursors. The model can be 
adapted to many categories of patients, or healthy platelet donors. It can also 
be modified to fit species other than human. By providing specific parameters 
one can adjust the model so as to yield particular predictions about the 
thrombopoietic profile of an individual patient. Platelet disorders, such as 
cyclic thrombocytopenia, may also be simulated. 



EXAMPLE 2: Neutrophil Bone Marrow and Peripheral Blood Compartment 
under the Effects of Growth-Factors and Treatment with Granulocyte 
Colony Stimulating Factor (G-CSF) 

The neutrophil lineage originates in pluripotent stem cells that 
proliferate and become committed to the neutrophil lineage. These cells then 
undergo gradual maturation accompanied with proliferation through the three 
morphologically distinguishable mitotic compartments: Myeloblasts, 
promyelocytes and myelocytes. The myelocytes then mature and lose their 
capacity to proliferate, and thus enter the post mitotic compartment. In the 
post-mitotic compartment the cells continue their gradual maturation, which is 
not accompanied with proliferation through the three morphologically 
distinguishable sub-compartments: Metamyelocyte, band and segmented- 
neutrophils. Cells exit the various sub-compartments in the post-mitotic 
compartment and enter the blood as neutrophils. They then migrate from the 
blood to the tissues. 

Th6 Granulocyte-Colony Stimulating Factor (G-CSF) effects an 
increase in blood neutrophil levels primarily by increasing production in the 
mitotic compartment and shortening the transit time of the post-mitotic 
compartment. 

Thus, the first compartment of the mitotic pool (myeloblast) receives an 
inflow of cells from stem-cell precursors. Inflow for each of the other 
compartments is from outflow of the previous one, subject to multiplication 
factors due to cell replication in the mitotic stages. 



Models regarding granulopoiesis in normal humans and in humans 
with pathologies of the bone marrow were suggested previously in order to 
give a coherent description of the kinetics of granulocytes from experimental 
data (Cartwright - 1964, Mary - 1978, Rubinow - 1974). In recent years 
Schmitz et al. developed a kinetic simulation model for the effects of G-CSF 
on granulopoiesis (Schmitz - 1993), and used it for the analysis of 
administration of G-CSF to patients suffering from cyclic neutropenia (Schmitz 
-1995). However, the data Schmitz rests upon for his model has been more 
accurately assessed in recent years by Price et al. (1998) and Chatta et al. 
(1996). Actual empirical data regarding compartment sizes and their transit 
times was not incorporated into their model despite the importance of these 
data (Dancey et al. 1976). 
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I. MODEL OF NEUTROPHIL LINEAGE AND EFFECTS OF G-CSF 
A. G-CSF 

The effects of G-CSF on the neutrophil lineage are relayed in the 
model in three stages. The first is the administered amount of cytokine given 

5 at time t, which is marked: G adm 

The Gadm vector serves as the control variable for optimization of G-CSF 
administration. 

The second stage represents the pharmacokinetic behavior of G-CSF 
in circulation. It incorporates, for instance, the half-life of G-CSF, and could in 
10 the future be modified to express more of the effects of time on G-CSF 



activity. This level is marked G[ 



blood ' 



G-CSF is eliminated from the blood in a Poissonic manner according to 
the following equation, as stated by Stute N, Furman WL, Schell M and Evans 
WE in "Pharmocokinetics of recombinant human granulocyte-macrophage 
15 colony stimulating factor in children after intravenous and subcutaneous 
administration"" Journal of Pharmaceutical Science, 84(7): 824-828, 1995: 



^ blood blood \ 1 7 >^ ^adm 

(14) 

where ty % is the half-life of G-CSF in the blood, and G\ lood = 

Recent data by Terashi K, Oka M, Ohdo S, Furudubo T, Ideda C, 
20 Fukuda M, Soda H, Higuchi S and Kohno S, in "Close association between 
clearance of recombinant human granulocyte colony stimulating factor (G- 
CSF) and G-CSF receptor on neutrophils in cancer patients", Antimicrobial 
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Agents and Chemotherapy, 43(1): 21-24, 1999, points to the dependence of 
the half-life of G-CSF on neutrophil counts. In the absence of exact kinetics 
of G-CSF effects on the neutrophil lineage, the half-life is considered as a 
constant, though this could be modified should more exact information 
emerge. 

Only exogenously produced G-CSF is considered to affect the kinetic 
parameters, and endogenously produced G-CSF levels and effects are set to 
zero. If more empirical data regarding the production of endogenous G-CSF 
is made available, it could be incorporated into the equation as well. 

The third and final stage models the pharmacodynamic effects of G- 
CSF on the kinetic parameters. As will be elaborated subsequently, the 
dependence of the various kinetic parameters of the neutrophil lineage on the 
level of G-CSF in the blood is assumed to be through either non-decreasing 
concave or non-increasing convex functions. This reproduces the effects of 
saturation that are seen in clinical studies on the effects of G-CSF, such as 
the study by Duhrsen U, Villeval JL, Boyd J Kannourakis G, Morstyn G and 
Metcalf D in "Effects of recombinant human granulocyte colony-stimulating 
factor on hematopoietic progenitor cells in cancer patients", Blood, 72(6): 
2074-2081, 1988. That is, addition of G-CSF carries a lesser effect when its 
level in circulation is already high. 

B. Biological Model 

Mitotic Compartment 

Long-term effects of G-CSF administration take place in the 
mitotic compartment. Although the major contributor to heightened blood 
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neutrophil counts in the short term is the post mitotic compartment's 
shortening of transit time due to G-CSF administration, this high level cannot 
be maintained over the long term without increased production in the mitotic 
compartment. 

The mitotic compartment is divided into subcompartments. The kth 
subcompartment contains all cells of chronological age between and k 
hours, relative to the time of entry into the mitotic compartment. The number 
of cells in subcompartment /rat time f is marked as m f k . 

10 

ke {L.t} 

(15) 

where t is the transit time of the entire mitotic compartment, and is assumed 
to be the same and constant for all cells entering the mitotic compartment, 
and U is a vector reflecting the flow of newly committed cells into the mitotic 
compartment. The biological grounds for this definition is the existence of a 
myeloid stem cell reservoir, which is known to supply new committed cells to 
the mitotic compartment. However, the reservoir's actual kinetics are not very 
well explored empirically. We therefore fix /, to levels such that the overall 
size of the mitotic compartment as well as the kinetics of the neutrophils in 
circulation would match those obtained empirically. 

Any new biological data that emerges may help define the kinetics 
more accurately within the framework of this model, although results of this 
model indicate that the assumption of a constant rate of stem cells flowing 
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into the mitotic compartment in the absence of G-CSF is plausible. For every 
ne{l..x}, and for every f, amplification occurs at the exit from m n \ according to 
Equation 1 5 as follows: 

<\\=<^ n (G t blood ) 

where: 

a n is a non-decreasing concave function of G-CSF levels in the blood, 
which determines the factor of amplification in the hourly subcompartment n. 
If, for instance, no amplification occurs at subcompartment n 0 at time X? then 

a = 1 
^^Gliood l ^( x n(G t blood )<2 

(17)The size of the morphological sub-compartments in the mitotic 
compartment at time t is determined as: 

(18)Where n 0 is the first hourly sub-compartment of a morphological sub- 
compartment and rii is its last hourly sub-compartment. The division into the 
morphological sub-compartments is used only for fine-tuning of the kinetic 
parameters with the use of experimental data. 

The mitotic compartment was modeled with an intention to facilitate the 
specific cell-cycle cytotoxic effects of chemotherapy. Therefore, cohorts of 
one hour are modeled as undergoing a process of maturation and 
amplification culminating in their entry into the post-mitotic as described 
below. Effects of chemotherapy may be incorporated into the model by 
mapping the various cell-cycle phases (G1, S, G2, M) to the hourly cohorts 



modeled and formulating a function of the cytotoxic effects of chemotherapy 
on these phases. 

The experimental literature shows wide agreement regarding the 
steady state normal amounts of circulating neutrohpils, size of the post-mitotic 
5 compartment and the three morphologically distinct sub-compartments of the 
mitotic compartment, and post-mitotic transit time and amplification rates in 
the mitotic sub-compartments (see, for example, Dancey JT, Deubelbeiss KA, 
Harker LA and Finch CA, in "Neutrophil kinetics" in Man. Journal of Clinical 
Investigation, 58(3): 705-715, 1976; Price TH, Chatta GS and Dale DC, 

10 "Effect of recombinant granulocytee colony-stimulating factor on neutrophil 
kinetics in normal young and elderly humans", Blood 88(1): 335-340, 1996; 
and Dresch Mary in "Growth fraction of myelocytes in normal human 
granulopoiesis", Cell Tissue Kinetics 19: 11-22, 1986). To determine other 
relevant kinetic parameters, which were either not available in the literature or 

15 were given a wide range by experimentalists, steady state kinetics was 
assumed and an iterative process was employed. These parameters include 
the inflow of stem cells to the myeloblast compartment and the transit times of 
the mitotic sub-compartments. 

The half life of blood neutrophils and the steady state number of 
20 neutrophils were taken as 7.6h and 0.4x1 0 9 cells/kg body weight, respectively 
(Dancey et al., 1976). Similarly, the same calculation may be made for each 
patient that is to be modeled. This would allow the dynamics of every patient 
to be described by the simulation. The average size of the post-mitotic 
compartment (5.84x1 0 9 cells/kg body weight - Dancey, 1976) and the transit 
25 time of the compartment (160h - Dancey, 1976; Mary, 1986; Price, 1996) are 



compatible with the size and half-life of the circulating neutrophil compartment 
reported by Dancey, thus supporting the steady state analysis. 

In order to determine the amount of cells in the hourly sub- 
compartments in the mitotic compartment, all compartments in the lineage 
5 were modeled using a steady state assumption. The number of cells exiting 
the circulating neutrophil pool equals the number of cells exiting the post 
mitotic compartment, which in turn equals the hourly production of cells in the 
mitotic compartment. Thus, the number of cells in the last hourly cohort of 
the mitotic compartment can be determined from the neutrophil decay rate, 

10 which is available in the literature. However, this calculation is based on 
assumptions that there is no apoptosis in the post-mitotic compartment. 
Direct experimental data by Thiele J, Zirbes TK, Lorenzen J, Kvasnicka HM, 
Scholz S, Erdmann A, Flucke U, Diehl V and Fischer R, in "Hematopoietic 
turnover index in reactive and neoplastic bone marrow lesions: Quantification 

15 by apoptosis and PCNA labeling," Annals of Hematology 75(1-2): 33-39, 
1997, suggests that apoptosis is not a significant phenomenon in normal 
human bone marrow. The size calculated for the mitotic compartment is 
close to that experimentally obtained by Dancey and Price, thus supporting 
the notion that apoptosis is not a significant kinetic factor in the lineage. 

20 Values for the production of cells in the mitotic compartment can later be 
modified in light of new evidence. 

Regarding the transit time of the mitotic compartment there is little 
agreement in the literature, with a range of 90-160 hours given by most 
experimentalists (see Dresch Mary in "Growth fraction of myelocytes in 
25 normal hman granulopoiesis," Cell Tissue Kinetics 19: 11-22, 1986). In order 



to determine the transit times of the mitotic morphological sub-compartments, 
as in Equation 18, the following constraints were considered: 

1. The sizes of the theoretically obtained morphological sub- 
compartments must fit those reported experimentally in normal 
human hematopoiesis (Dancey, 1976) and under the effects of G- 
CSF (Price, 1996); 

2. At least 24 hours, the typical cell cycle, must separate 
amplification points; 

3. The size of the last hourly sub-compartment must equal the hourly 
production of the mitotic compartment (calculated with the 
aforementioned iterative process assuming steady state kinetics); 

4. Amplification inside the compartment is set at the levels 
determined by Mary 1 984; and 

5. The total transit time of the mitotic compartment must be within 
the 90-160 hour range. 

By using the values shown in Table x, an excellent fit was obtained within the 
above-mentioned constraints. 

It should be noted that when other alternatives with shorter transit 
times were attempted, results could not be obtained that agreed with the 
literature regarding the size of the mitotic pool or its production. 
Furthermore, a fit between our simulation model's results regarding 
PMN counts in peripheral blood with empirical data could not be 
achieved without speculating extensively on the nature of G-CSF effects 



on non-committed stem cells (data not provided). It should be noted, that 
little empirical quantitative data is available regarding stem cells. 

The effects of G-CSF on this compartment are modeled as an increase 
in the rate of cells entering the myeloblasts from the uncommitted stem cell 
5 pool, increases in the rates of mitosis, and introduction of new points of 
amplification as shown in Equation 15, 16. Since little data is available 
regarding the increases in amplification due to G-CSF, an initial assumption 
was made that amplification reaches full potential at points that under normal 
conditions undergo an amplification of below a factor of 2. Additionally, it was 
10 assumed that the transit time in all mitotic sub-compartments and the typical 
cell cycle duration are not affected by G-CSF, based on lack of evidence to 
the contrary. 

Reference is now made to Figs. X and y, which show a comparison of 
neutrophil production according to the described model and experimental 
15 data in the literature. Increased neutrophil production is in accordance with 
the neutrophil counts reported by Price et al. In addition, these increases are 
in accordance with Price's data about neutrophil bone marrow pool sizes. 

Reproduction of the effect of G-CSF on neutrophil counts and the 
mitotic compartment sizes beyond day 5 of administration was accomplished 

20 by assuming an increase (15% with the highest dose of G-CSF) in the rate of 
cells entering the myeloblast compartment. Alternatively, G-CSF may change 
the behavior of the myeloblast compartment such that some of the cells there 
undergo self-renewal instead of moving on to the promyelocyte compartment. 
However, no empirical data to support this is available. The model can be 

25 modified in light of new experimental data in the future. 
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Post-mitotic compartment 

The post mitotic compartment is relatively insensitive to cytotoxic 
chemotherapy. Therefore, it is biologically acceptable and computationally 
sensible to model this compartment as a single pool of cells, such that the last 
hourly cohort of the mitotic compartment enters the compartment, and a 
proportion of the cells within the compartment enter the neutrophil pool every 
hour. 

The post mitotic compartment at time t is a single quantity of cells p\ 
such that: 

(19) 

where l 3 is a convex, non-increasing function of G-CSF levels in the blood, 
which takes values in the range of [0-1]. This definition entails p'>0. 

G-CSF affects the post-mitotic compartment by shortening its transit 
time (i. e. decreasing l 3 ). Price notes that the number of cells in the post 
mitotic compartment is not significantly changed following administration of G- 
CSF. This determination is based on counts made on day 5 after G-CSF 
administration. Thus, it can be safely assumed that any increased production 
of the mitotic compartment flowing into the post mitotic compartment is 
translated over the long-term to an increase in the flow of cells from the post 
mitotic compartment to the neutrophil pool. This increased flow is 
compensated by increased production in the mitotic compartment only at a 



later stage. Therefore, an upper limit to the number of cells in the post mitotic 
compartment was set, which is at the values given as steady state counts (n). 

In brief, the effects of G-CSF on the neutrophil lineage are modeled 
during the first few days primarily as a decrease in the counts of the post- 
5 mitotic compartment, which is then compensated by an increased production 
in the mitotic pool. This compensation sustains the increase in neutrophil 
counts in peripheral blood. 

Reference is now made to Fig. 1 1 , which is a graphical illustration of a 
simulation of the model. Though no empirical data is available on this point, 

10 simulations of the model predict that the number of cells in the post-mitotic 
compartment decreases substantially during the first two days of G-CSF 
administration, and then replenishes, so that on the sixth day the counts 
return almost to their normal levels. This replenishment lags behind that of 
Price et al report by a few hours. We can thus formulate a testable 

15 hypothesis, i.e., whether using the same G-CSF protocol Price et al used, 
there is indeed a nadir on day 3 of the treatment. 



Compartment 


Day 0 (no G- 
CSF) (xlO 9 
cells/ kg. 
Body 
weight) 


Day 15 of G-CSF 
treatment (xlO 9 cells/ 
kg. body weight) 


Relative increase in 
compartment size 
due to G-CSF 


Myeloblasts 


0.140 


0.153 


1.09 


Promyelocytes 


0.582 


0.898 


1.54 


Myelocytes 


1.373 


3.564 


2.60 


Mitotic Total = 


2.10 


4.615 


2.20 


Circulating neutrophils 


0.4 


2.35 


5.88 



Table 1. Simulated kinetics after 15 days of subcutaneous administration of 300ug G-CSF/kg weight. Day 0 
values are the mean values Dancey et al (1976) use. 



20 We shall mark as o* the outflow from the post mitotic compartment: 

t t t ?+i 
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(20) 

The number of neutrophils in the circulating blood compartment at time t is 
marked n* and is modeled as a single quantity of cells, such that: 



n t+l =o f +wYl-— ) 

t ^ 2 (21) 

where t Vz is the half-life of neutrophils in the blood, as defined in the biological 
literature, and is assumed to be constant regardless of G-CSF levels (Lord 
1989), though this could be easily modified. The kinetics of neutrophils in the 
tissues are not modeled in this work. 



Neutrophils and G-CSF in the Circulating Blood 

The elimination of neutrophils from peripheral blood follows a Poisson 
distribution, and can therefore be described as an exponential function 
(Cartwright, 1964). Therefore the rate of cells leaving this compartment is 
is based on half-life determinations available in the literature. Since no direct 
cytotoxic effects of chemotherapy have been described for this compartment 
it is also modeled as a single pool of cells. 

At the normal healthy level we have the following relationship: 



& = m. '=()' =-^xln2 
T 

which reflects the stability of the steady state. 



The kinetics of G-CSF is also modeled as an exponential distribution 
with a half-life of 3.5 hours (Stute, 1992) (Eq. 14). 

The effects of G-CSF on the kinetics of the neutrophil lineage appear 
not to be a linear function of G-CSF administration levels. Since data provided 
in the literature (Chatta 1994) only refers to two doses (30 and 300^igram/kg 
body weight) we can only speculate on the effects of other levels of G-CSF. 
After trial and error analysis, it was found that assuming that the effects of the 
300-y.gram dose are the maximal, at the 30-^igram its effects are about 30% 
of the maximum. 

Reference is now made to Figs. 12A and 12B, which are graphical 
illustrations of the effects of G-CSF at the two doses. The effects as a 
function of G-CSF level are connected piece-wise linearly. This way, the 
neutrophil levels observed clinically under both the 300 and the 30-fxgram 
protocols are obtained. 



II. LINEAR IMPLEMENTATION OF THE MODEL 

This model will later be incorporated into an optimization scheme that 
will have as its objective function both the aims of minimizing G-CSF 
administration and returning the neutrophil lineage to its normal levels. 

Although the above-outlined model may be implemented in any 
number of optimization methods, linear programming was chosen because of 
its inherent advantages compared with other techniques, i.e. its ability to 
provide an optimal solution using partially analytical methods, and therefore 



being more computationally tractable (Gill 1991). On the other hand, 
implementation of this model in linear programming carries with it the 
disadvantage that certain computations must be approximated linearly since 
they cannot be performed directly using linear methods. Thus, we shall 
5 compare the 'closeness' of the solution obtained through linear programming 
will be compared with that obtained through another, non-linear method of 
optimization. 

The significant parts of the model that must be modified due to the 
linear programming implementation are the sections in which multiplication of 
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io two variables is defined, since this operator is not itself linear. Therefore, 
multiplication is defined as an approximated value constrained within 
piecewise linear constraints that most closely bound the product within a four- 
faced polyhedron in 3-dimensional space whose vertices are 



is Where x mint x max , y m \n, Vmax are the constant biologically defined minima and 
maxima of x and y. 
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(22) 



Multiplication may also be approximated with variations on the linear 
least squares method, by finding one plane that's closest to the four vertices 



20 defined. 



The other functions that need to be defined linearly are those 
concerning the pharmacodynamics of G-CSF. Due to the nature of these 
functions (either non-increasing convex or non-decreasing concave), these 
effects are implemented as piece-wise linear functions whose breakpoints are 
the doses for which actual experimental data is available (Chatta 1994). Note 
that the effects of G-CSF on each of the kinetic parameters have not been 
determined in a detailed manner by experimentalists. Rather its effects over 
a few dose levels on the neutrophil blood counts and the size of the 
morphologically different mitotic compartments and the post mitotic 
compartment have been determined. From these data, the effects of G-CSF 
on the actual kinetic parameters (probability of mitosis, transit time and inflow 
of cells into the myeloblast compartment from stem cell progenitors) has been 
reconstructed at the dose levels available in the literature. These points are 
then connected linearly to obtain piecewise linear functions relating G-CSF 
levels to their effect on those parameters. Further experimental data in the 
future could be used to produce more accurate functions. 

At the amplification points within the mitotic compartment, the linearly 
approximated multiplication operator (Eq 22) is used instead of the product 
defined in Eq. 16. 

At points where no amplification occurs the quantity from one 
compartment is simply transferred to the next according to the following 
Equation: 



(23) 



Values are set according to the steady state values of the mitotic 
compartment, or are depleted according to the kill function of the 
chemotherapy. 



defined as a linear approximation of a product. 

3. Formulation of the Model as an Optimization Problem for Linear 
Programming 

The simulation spans a finite number of discrete time steps denoted by T. 

10 We define as the control variable the vector that represents G-CSF 
administration at every given laous t Tl 



The objective function is defined as maximization of the following expression: 



is where p f is the number of cells in the post mitotic compartment at time t, and 
• is a scalar weighting coefficient. The logic for formulating the objective 
function this way is that the ability to maintain the post mitotic compartment's 
steady state size for a prolonged period of time is sufficient for rehabilitation 
of the neutrophil lineage as a whole. Also, our goal is to minimize the total 

20 administered quantity of G-CSF. • is introduced to allow us to factor in both 
these goals in one objective. Also, this would allow a different weight to be 
given for certain times, e.g. were it determined (by clinicians) that the later 
states of the post-mitotic compartment should be weighted more than the first 
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The flow out of the post mitotic compartment (Eq. 20) is similarly 



T 




(Eq. 24) 



ones. Obviously this is only one of the possible formulations of the objective 
function as defined in the previous section. 

The pharmacokinetics and pharmacodynamics of G-CSF that were 
defined generally in the mathematical model are defined piecewise linearly. 
Some of the considerations that we put into formulating these functions were 
based directly on experimental evidence (elaborated in the main body of text). 
We note however, that actual experimental data regarding the direct effects of 
G-CSF on the kinetic parameters in which this model is interested is rather 
scant. Therefore, some formulations were conducted through partly analytic 
and partly trial-and-error methods. 

The formulation of the model in piecewise linear terms will allow use of 
this model as a clinical tool in three ways. First, the model will determine the 
effectiveness of various protocols suggested by clinicians prior to their actual 
use on human patients. Second, the model allows computation of the optimal 
protocol in a given situation of neutrophil counts, so that the neutropenic 
period following chemotherapy is either shortened or completely avoided at a 
minimal cost and exposure to G-CSF. Third, the model serves as a 
constituent in a broader framework of clinical tools that will compute the most 
optimal treatment plan for chemotherapy and growth factors. These uses 
should help clinicians administer more rational treatment to their patients 
minimizing both suffering and medical costs. 



Amplification at 


Mean transit time 


Size (10 6 cells/kg 


Compartment 


the exit 


(hours) 


weight) 




2 + 


24 _ 


0.139* 


Myeloblasts 


2 + 


48" 


0.558* 


Promyelocytes 


1.5 + 


48- 


* 

1.4 


Myelocytes 


1 


160* 


5.84* 


Post mitotic 


0 


10.96* 


0.4* 


Neutrophils 



Table 1. Kinetics under steady state conditions in healthy humans. 

*Dancey, 1976. 
+ Dresch, 1986. 

"Calculated based on the steady state assumption as elaborated in the main text. 



EXAMPLE 3: Cancer and Treatment with Cytotoxic Drug Delivery 

Introduction 

Cancer is the second leading cause of mortality in the US, resulting in 
approximately 550,000 deaths a year. There has been a significant overall rise 
5 in cancer cases in recent years, attributable to the aging of the population. 
Another contributing factor to the rise in the verifiable number of cases is the 
wider use of screening tests, such as mammography and elevated levels of 
prostate specific antigen (PSA) in the blood. 

Neither better detection nor the natural phenomenon of aging, 
10 however, can entirely explain the increase in new cases of tumors. Meanwhile, 
other cancers, like brain tumors and non-Hodgkin's lymphoma, are becoming 
more common. Their increase could reflect changes in exposures to as yet 
unidentified carcinogens. Current trends suggest that cancer may overtake 
heart disease as the nation's no. 1 killer in the foreseeable future. As gene 
is therapy still faces significant hurdles before it becomes an established 
therapeutic strategy, present control of cancer depends entirely on 
chemotherapeutic methods. 

Chemotherapy is treatment with drugs to destroy cancer cells. There 
are more than 50 drugs that are now used to delay or stop the growth of 
20 cancer. More than a dozen cancers that formerly were fatal are now treatable, 
prolonging patients' lives with chemotherapy. 

Treatment is performed using agents that are widely non-cancer- 
specific, killing cells that have a high proliferation rate. Therefore, in addition to 
the malignant cells, most chemotherapeutic agents also cause severe side- 
25 effects because of the damage inflicted on normal body cells. Many patients 



develop severe nausea and vomiting, become very tired, and lose their hair 
temporarily. Special drugs are given to alleviate some of these symptoms, 
particularly the nausea and vomiting. Chemotherapeutic drugs are usually 
given in combination with one another or in a particular sequence for a relatively 
5 short time. 

Chemotherapy is a problem involving many interactive nonlinear 
processes which operate on different organizational levels of the biological 
system. It usually involves genomic dynamics, namely, point mutations, gene 
amplification or other changes on the genomic level, which may result in 

10 increasing virulence of the neoplasia, or in the emergence of drug resistance. 
Chemotherapy may affect many events on the cellular level, such as cell-cycle 
arrest at different checkpoints, cell transition in and out of the proliferation cycle, 
etc. Chemotherapy may also interfere with the function of entire organs, most 
notably, with bone marrow blood production. In recent years molecular biology 

15 and genetics has made an important step forward in documenting many of 
these processes. Yet, for assessing the contribution of specific molecular 
elements to the great variety of disease profiles, experimental biology must be 
provided with tools that allow a formal and systematic analysis of the intricate 
interaction between the genomic, cellular and cell population processes in the 

20 host and in the disease agent. This system is so complex that there is no 
intuitive way to know how small changes in the drug protocol will affect 
prognosis. But in spite of this intricacy, attempts to improve chemotherapy have 
been carried out by "trial and error" alone, with no formal theory underlying the 
application of specific drug schedules. Such an approach "is apt to result in no 



improvement, only discouragement and little useful information for future 
planning" (Skipper, 1986). 

The treatment of cancer by cytotoxic drug (or drug combination) 

5 delivery is addressed. In this model, two generic types of cells are considered: 
host cells and target cells. Target cells are, in fact, the tumor. Both types of 
cells may be damaged while exposed to chemotherapy. The aim is to obtain 
the most suitable treatment protocol according to specified conditions and 
limitations. It is assumed that the cell dynamics are deterministic and known, 

10 and that both types of cells are sensitive to chemotherapeutic agents in certain 
known period of the cell-cycle phases. If a cell is exposed to chemotherapy 
during part of its critical phase, there is a chance that it will be eliminated, 
blocked or affected in any other known way. The description of the dynamics of 
the delivered drugs is assumed to be known as well. 

15 In order to achieve the goal optimization process is applied to the 

model. The optimization module uses the model predictions in order to search 
for the suitable solution to posted optimization problem. Precise defining of 
optimization objectives as well as the relevant parameters adjustment is done 
according to the settings defined by user/operator for every special case. The 

20 method can be applied in general cases as well as in specific ones. 

2. Model of Biological System 

The basic layer of the model incorporates a description of age 
distribution of cycling cells and number of resting (quiescent) cells. The term 



"age of the cell" here refers to chronological age starting from the 
conventional beginning point of mitotic cycle. 

Reference is now made to Fig. 1 3, which is a schematic illustration of 
the tumor cell cycle layer. The whole cycle is divided into 4 compartments, or 

5 stages (Gi, S, G 2 and M). Each compartment is divided into equal 
subcompartments, where i th subcompartment in each stage represents cells 
of age i in the particular stage (i.e. they have spent i time-steps in this stage). 
The quiescent stage is denoted G 0 . The cell cycle follows a direction as 
shown by arrows (#). Thus, cells enter each stage starting with the first 

10 subcompartment, denoted Gi . 

The model can be described mathematically as follows: Let T k denote 
the maximum duration of k^ stage in the cycle. Let • t denote the time 
resolution of the model in discrete time steps. X' k (t) is a function, which 
represents the number of cells in stage k in the I th sub-compartment, at time f 

15 to f+» t. Both time and age are measured in the same unit, in this case, 
hours. Let Q(t) represent the number of resting cells at time r to t+» t. 
Trans(k,i,t) represents the probability that a cell of age / in the stage k will 
move to the next (k+1) compartment. Cells entering the new stage always 
start from the first subcompartment, i.e. from /=1 . This probability may change 

20 with time, representing the influence of conditions on cycle length distribution. 

By definition, the cell cannot remain more than T k timesteps in the \c 
compartment, as described in the following equation. 

V k : Trans ( k , T k ) = 1 
The restriction point (R-point) represents a cell's commitment to 

25 complete the mitotic cycle. Let T R denote the age at which the cell passes 



through the restriction point in G1. Only cells in Gi with l<T R can the cycle 

(in the absence of a drug). 

The total number of proliferating cells P(t) can be calculated as follows: 



k = G 1 ,S ,G 2 ,M 



( T„ 

i = 1 



V 



j 



In every time interval, quiescent cells may return to the proliferation 
pool. Alternatively, proliferating cells may change their state to become 

10 quiescent if and only if they are in the G1 stage and at age /, where T R * i>0.. 
To describe this process we introduce the function Gi. 0 (i,t) which describes 
the number of G1 cells in age / which become quiescent during time interval 
[t, t+ • fj. This function may receive negative values, accounting for cells that 
return from resting to proliferation. 

15 As we assume that the exit to quiescence can occur only prior to the R- 

point (even in cancer cells), and that a resting cell that returns to proliferation 

i 

enters the cycle at T 0 , it can be stated: 

Vi>T R ,Vt:G^ 0 (i,t) = 0 

20 It must be noted that this function is not dependent on /' and r solely. Its value 
is determined according to current cell distribution and all the general 
parameters that characterize the described cells group. The same should be 
said about the values of Trans vector that can change during the history of 
given population. 



The model traces the development of described group of cancer cells 
using given parameters, by calculating the number of cells in each and every 
subcompartment according to the following stepwise equations: 




10 



15 for k=G 2 ,M, k-1 returns the previous stage (e.g. G 2 -1=S). 

These equations make it possible to calculate the number of cells in 
each subcompartment at every time interval [t* t] starting from initial 
distribution (e.g. at time t=0). Since in this model cell ages are measured in 

20 absolute time units, these measurements refer to the chronological age of the 
cell, and not the biological age, whose units are relative to a maturation rate 
that differs from cell to cell. Consequently, in this model no cell can remain in 
the same age subcompartment after every time step. On the other hand, a 
fraction of the cells that leaves any subcompartment may be transferred to 

25 the first subcompartment of the next stage, according to probability vector 

sn 



Trans(k,i,t). This vector provides the ability to account for variability of cycle 
lengths while retaining a deterministic approach. 

The behavior of the cell population in this model is completely 
controlled by two components: Trans vector, and G1 ->G0 function. These two 

5 functions determine uniquely the outcome of every single time step, and, 
consequently the result over long periods. Thus, they are referred to as 
"control function^'. The values of these functions may be dependent not only 
on time and age of cells, but also on the current population state (or, 
generally, on the whole history of the population) as well as on the 

10 environment associated with a given cell group. However, those parameters 
are similar for all the cells in the group, implying that the model presented 
here is suitable for describing highly homogenous group of cells. Therefore, 
the basic layer of the model should give a realistic description for a uniform 
group of cancer cells for which environmental conditions and relevant 

15 biological properties are defined, in a way that will allow the construction of 
the control functions for the group. 

2.2. The General Tumor Model 

20 In the general approach the whole model is viewed as constructed from 
similar components, each of them derived from the basic structure described 
in the previous paragraph. Each component represents cells that are 
subjected to the same environmental conditions and, therefore behaves 
similarly (to be denoted homogenous group). The whole tumor is modeled 

25 as a union of many varying homogenous groups of cells, where the 



development of each group can be accurately predicted (when local 

conditions are known). 

This general model simulates progress of the tumor in discrete steps of 
time. At each step the number of cells in each subcompartment of each 
group is calculated according to the previous state, parameters of tumor, drug 
concentration, etc. The parameters of the tumor must include all the 
information that is relevant to prognosis. Some of these parameters are 
defined locally, e.g., those relating to the tumor's geometry. For this reason 
the representation of the spatial structure will be included. 

The cells will be able to pass between the groups during the 
development of the tumor. This allows the representation of the changes in 
the local conditions during the tumor evolution (e.g. forming of necrotic core, 
improvement in "living conditions" in vascularized regions, etc.). In addition, 
all the parameters of the tumor may change in accordance with the dynamics 
of the cancer. 

The calculation of the tumor development over time will be done by 
stepwise execution of the described simulation and can be used to predict the 
outcome of the treatment or in fitness function for search algorithms. 

2.3. From General to Individual Tumor Model 

When the general theoretical description of the model is accomplished, 
the model is fitted to represent the actual tumors. We render it patient-specific 
by adjusting all the parameters that determine the behavior of the modeled 



tumor to those of the real cancer in the patient's body. In order to accomplish 
this task we will establish the connections between mathematical parameters 
(most of them will have direct biological implication) and every kind of data 
that is practically obtainable in the clinic. These connections may be defined 
5 through research on statistic correlation between different parameters 
(including genotype-phenotype correlation), or using advanced biochemical 
research (which may establish quantitative relations). 

Thus defined, the model will be able to give realistic predictions for 
treatment outcomes either for specific patient or for a broad range profiles of 
10 patients and diseases. This tool can serve to perform the prognosis of either 
an untreated cancer patient, or as a basis for treatment modeling as is 
described below. 

15 3. Introducing Pharmacology 

In order to simulate cancer treatment we add to the above model the 
pharmacologic component. We model pharmacokinetics as well as 
pharmacodynamics for specific anticancer drugs. We begin by representing 
20 cell-cycle specific drugs, however the model implies no restriction on the type 
of drugs to be employed. 

We model the distribution of the drug in and around the tumor as well 
as in the blood (the drug kinetics). For this purpose, we use the suitable 
model, defining it precisely for every certain type of the drug. The 
25 concentrations of drug in the body are calculated at every time step in 



accordance with the drug administration specified by the protocol, and 
different processes that define drug kinetics in the body. 

The dynamics of the drug are represented through the direct influence 
of the drug on tumor cells. The effects on the proliferating cells are mostly 
blocking the cycle in different stages (which can be modeled as cell arrest) 
and cell death (immediate or after being in the block). Cell-cycle specific 
drugs are believed to have no direct influence on quiescent cells, but can 
affect them indirectly by killing proliferative cells and therefore changing local 
conditions. Where additional types of drugs added to the model, their effect 
on any kind of cells may be too modeled as killing certain fraction of cells 
(which is dose-dependent) or changing the behavior of the cells. 

Additional phenomena that may prove significant in drug kinetics and 
dynamics (e.g. rate of absorption by the cells, development of tumor 
resistance to specific drug, etc.) can be introduced into the model to make it 
as realistic as needed. 

The description of the drug in the model is done in terms of quantitative 
functions, which enable to calculate the drug amounts at certain locations and 
the tumor response to it at every time step. In the general case, these 
functions include parameters that depend on the specific data (drug type, 
body parameters, characteristic of the tumor, etc.) and can be determined in 
given situation (patient/case). The relation between clinical data and these 
parameters can be established in the ways similar to those described for the 
cancer model. 

The combination of cancer model with the drug model described above 
makes it possible to predict the outcome of the treatment, given the relevant 



parameters for the drug, the cancer and the patient. Again, the prognosis may 
be made for specified cases as well as for broad profiles of patients or 
disease. This simulation also serves to build the fitness function used for the 
optimaization objectives. 

5 

4. Combining with minimizing host toxicity. 

Although an accurate predictive tool, the model that represents 
10 chemotherapy of tumor alone cannot be used in optimization, for it posts no 
constraint on the choosing of the treatment. Actually, this model implies using 
as much drug as possible until the final elimination of the tumor; while in the 
live system the toxicity of the drug is the most important constraint limiting the 
treatment. In most cases of anticancer chemotherapy the dose-limiting 
is toxicity is bone marrow suppression, the two most sensitive bone marrow 
lineages being granulopoiesis and thrombopoiesis. Accordingly, those two 
were chosen as an example and are modeled separately and in a similar way 
to predict the negative effect of the chemotherapy on them. These models 
reconstruct the damage caused by the chemotherapy to the bone marrow 
20 cells and the recovery of these lineages (treated by specific growth factors). 

Thus, the whole system is capable of predicting the result of 
chemotherapy treatment for the tumor as well as for bone marrow cells, 
allowing the use of the protocols that combine anticancer drugs and growth 
factors for healthy cells. 

6^ 



Chemotherapy toxicity to any other normal host cell population can be 
similarly taken into account, if it is defined as relevant for dose and schedule 
optimization. 

5 

VI. INDIVIDUALIZATION OF THE MODELS 

Due to a great degree of heterogeneity between malignant tumors (even 
among similar tissue types) and between patients, it would be advantageous 
to adjust the treatment protocol to the individual case. This individualization 
10 procedure includes three aspects: 

1 ) individual parameters of tumor dynamics 

2) individual parameters of patient-specific drug pharmacokinetics 

3) individual parameters of the dynamics of dose-limiting normal host 
tissues. 

15 Relevant data concerning individual cases can be obtained by research 

on statistic correlation between different parameters (including genotype- 
phenotype correlation), or using advanced biochemical research (which may 
establish quantitative relations). In the general model, important dynamic 
parameters are estimated from experimental studies conducted in certain 

20 patient populations. Any of these parameters, when available on the per 
patient basis, can be individualized, while those that are unavailable can be 
left as a population-based figure. This approach allows continuous increase 
in the degree of individualization of the treatment protocols with progress in 
the technology of parameter evaluation (e.g., oncochips). 

6? 



All different parameters may then be adjusted, which will result in an 
adjusted array of models to be simulated. Parameters may include many 
different factors, which are adjustable according to the needs of a 
pharmaceutical company for general use of the treatment, or may even be 
5 individualized for use by a specific clinician for a particular individual. Examples 
of parameters may include, but are not limited to age, weight, gender, previous 
reaction to treatment, desired percentage of healthy body cells, desired length 
of treatment protocol, pathologic or cytologic specifics, molecular markers, 
genetic markers etc. 

10 In order for the system to be user-friendly, all possible parameters are 

termed in ways that are easy for the user to understand. 

Once all the parameters are set, an array of solutions is produced 
based on the input parameters. A number of possible protocols can be set (is 
thus generated by the computer), a fitness function is applied, which results in 

15 scores for each of the proposed solutions. 

GENERATION OF PROTOCOL SPACE 

Refer back to Fig. 2. This model makes it possible to check any given 
20 treatment protocol and to choose a very good one according to user's criteria. 
The user may be a physician, a drug developer, a scientist, or anyone else 
who may need to determine a treatment protocol for a drug. The specific 
parameters may include several categories, such as individual patient 
characteristics and/or medical history, needs of a specific user (research, 



efficacy, treatment, etc.), and other particulars (such as maximum length of 
treatment, confidence level, etc.). 

. That is, an array of possible treatment protocols is created from which 
5 the optimal treatment protocol can then be chosen. It should be noted that the 
method does not imply the fitness estimation for all possible protocols. The 
use of operation research allows a much more sophisticated, yet resource 
saving procedure. 

An example of this procedure will be described as it relates to cancer 
10 treatment by chemotherapy, as described in the third embodiment of the 
invention above. However, it should be noted that a similar procedure may be 
performed in any of the embodiments. 

The procedure implements cell growth and cell death procedures, as 
defined in the detailed model above. There are certain pre-defined 
15 parameters, including the lengths of the host and target life-cycles, the 
lengths of their critical phases, and a resolution factor, that determines the 
length of a single time unit. The user is asked to define an action (treatment 
or non-treatment). Simulation of cell growth and death is then performed for 
the single time unit. This procedure is repeated until the end of the total 
20 simulation. Alternatively, the choice of treatment or non-treatment is made by 
the processor, with many possible permutations considered. In that case, the 
protocol space would be very large, and the resolution would depend on the 
(selected) length of the time unit (and computer capacity) 



fr5 



There are two procedures: one for growth simulation during treatment 
and one without. The array in which the numbers of cells are kept is updated 
once per time unit, whether with or without treatment present at that time. 

VII, Defining the Fitness Function 

The, fitness function is an important tool in Operation Research. In this 
case of protocol optimization, it allows the comparison between a number of 
different protocols each one of them scoring differently with respect to various 
objectives that can be set by the developers or by the users and identifying 
the protocol for which the highest weighted score is predicted. The fitness 
function calculates for any given protocol its relative efficacy ("score" or 
"fitness"), thus enabling a definite decision of the best protocol from a given 
set of protocols. 

In different cases, different objectives can be formulated. There are 
several settings in which such a model can be used, including but not limited 
to: 

One) clinical practice- where objectives can change depending on 

type of disease, condition of the patient, purpose of treatment, etc. 
Two) pharmaceutical company- where objectives can be aimed at 

finding the therapeutic window and an optimal schedule. 
Three) scientific setting- research oriented objectives can be aimed at. 
In each case, a particular fitness function can be formulated, reflecting 
all given requirements. Thus, in any particular case one can compare 



between different protocols and obtain the most suitable to his/her special 
purposes and needs. 

Examples for some alternative objectives are given: 

1 . The smallest number of drug dosings required for achieving any 
5 given aim. 

2. The lowest total drug dose required for achieving any given aim. 

3. The minimal total amount of drug needed for rehabilitation. 

4. The smallest deviation from the baseline at normal cell 
population count (e.g., platlet nadir) after chemotherapy or another cell- 

10 suppressive treatment. 

5. The shortest period of disease (e. g. thrombocytopenia). 

Using the fitness function it is possible to a) estimate the efficacy of a 
given protocol, b) search for the solution of an optimization problem, i.e., 
predict which protocol will be best of many potential protocols considered for 
15 curing/relieving the patient. 



VIII, Solving the Optimization Problem 

The optimization problem is stated using the described models: to find 
20 the protocol for drug administration (with option to growth factor 
administration) which maximizes the given fitness function. 

As explained above, we define the fitness function according to the 
user requirements. For example, the goal of the treatment may be defined as 
minimizing the number of cancer cells at the end of the treatment, minimizing 



the damage to the BM cells throughout the treatment or at its end, and curing 
the patient (where cure is defined precisely) as quickly as possible. Note that 
the fitness function may also include goals such as maximizing life 
expectancy, minimizing cost of treatment, minimizing treatment hazards 

5 and/or discomfort etc. Generally, the aim of optimization is to find the best 
protocol, i.e. the protocol that generates the best value of fitness. 
Customarily, this is achieved by mathematical analysis. However, 
mathematical analysis is restricted to over simplified models, whereas, in 
order to accommodate biologically realistic parameters, the described models 

10 are very complex and, therefore, cannot be solved analytically. On the other 
hand, the practical purpose of the treatment is not to find the best possible 
protocol (i.e., the global optimum) but "only" one that will suit the user's 
objectives, even if its fitness is not absolutely the best (i.e., the local 
optimum). For this reason we can be satisfied with the solution that can be 

15 shown to promise the pre-defined objectives. 

Hence, the optimization problem may be reformulated as follows: for 
given initial conditions, find the treatment protocol which will fulfill the user's 
requirements (e.g. curing a patient according to given definitions of cure) and 
subjected to given limitations (e.g. treatment duration, drug amounts, etc.). To 

20 this end it is not compulsory to find the global solution. It is enough, with 
regards to objectives and limitations, to perform search, using search 
algorithms, in certain regions of the protocols' space, and find the local 
maxima of the fitness function. After determining the locally best protocols, we 
can verify that they serve one's objectives and check them numerically for 

25 stability. 



Such a strategy can be used for identifying patient-specific treatment, 
as well as in the general case, where only the profile of the disease and the 
drug are specified. If more patient-specific data are supplied, the solution will 
be tailored more specifically. On the other hand, the optimization program 

5 could propose general recommendations for the protocol types for certain 
kinds of disease, treated with a certain kind of medication. 

It will be appreciated that the present invention is not limited by what 
has been described hereinabove and that numerous modifications, all of which 
fall within the scope of the present invention, exist. For example, while the 

10 present invention has been described with respect to certain specific cell 
lineages, the concept can be extended to any other lineage and treatment 
protocol which can be detailed mathematically (e.g., viral or bacterial diseases). 
Furthermore, certain assumptions were necessarily used in computing the 
mathematical models of the embodiments. Values and equations based on 

15 these assumptions can be changed if new information becomes available. 

It will be appreciated by persons skilled in the art that the present 
invention is not limited by what has been particularly shown and described 
herein above. Rather the scope of the invention is defined by the claims which 
follow: 



What is claimed is: [independent claims listed now, dependent claims to 
follow] (Note that we are not in a position to revise the claims. Z. Aaurt 

1 . A system for individualized optimization of a treatment protocol, the system 
comprising: 

a system model comprising: 

a model of a biological process; and 

a mathematical model of treatment effects on said biological 

process; 

a plurality of treatment protocols; 

means for adding individualized parameters to said system model, 

whereby said system model is modified based on said individualized 
parameters; and 

a selector for selecting an optimal treatment protocol from said 
plurality of treatment protocols based on modification of said system model. 

2. A system for identification of optimal treatment protocols using heuristic 
analysis, the system comprising: 

a biological model; 

a plurality of treatment protocols included within said model, thereby 
providing a plurality of treatment models, wherein said treatment models 
provide effects of said treatment protocols on said biological model; and 

heuristic means for identification of said optimal treatment protocol 
from said treatment models. 



What is claimed is: [independent claims listed now, dependent claims to 
follow] (Note that we are not in a position to revise the claims. Z. Aqur) 

1 . A system for individualized optimization of a treatment protocol, the system 
comprising: 

a system model comprising: 

a model of a biological process; and 

a mathematical model of treatment effects on said biological 

process; 

a plurality of treatment protocols; 

means for adding individualized parameters to said system model, 

whereby said system model is modified based on said individualized 
parameters; and 

a selector for selecting an optimal treatment protocol from said 
plurality of treatment protocols based on modification of said system model. 

2. A system for identification of optimal treatment protocols using heuristic 
analysis, the system comprising: 

a biological model; 

a plurality of treatment protocols included within said model, thereby 
providing a plurality of treatment models, wherein said treatment models 
provide effects of said treatment protocols on said biological model; and 

heuristic means for identification of said optimal treatment protocol 
from said treatment models. 



3. A system for presenting a individualized model of a biological process, the 
system comprising: 

a processor for generating a model of a biological process; 

a plurality of experimental data, included within said model; and 

means for including individual data into said processor, thereby 
providing said individualized model of said biological process. 

4. A system for prediction of outcomes of a treatment protocol, the system 
comprising: 

5. A system for describing treatment effects on a biological process 

6. A method for individualized optimization of treatment, the method 
comprising the steps of: 

providing a biological model; 

inputting parameters into said model; 

computing treatment protocols for said model within said parameters; 

applying a fitness function to said treatment protocols; and 

choosing an optimal treatment protocol based on said fitness function. 

7. A method for optimization of a treatment protocol, the method comprising 
the steps of: 

providing a variable representing an amount of medication 
administered; 

providing behavior characteristics of said administered medication; 



providing a model of a cell line predicting effects of said administered 
medication on said cell line; 

calculating a 

determining a treatment protocol 

5 8. * A system for prediction of outcomes of treatment of thrombocytopenia by 
Thrombopoietin (TPO), the system comprising: 

9. * A system for prediction of outcomes of treatment of neutropenia with 
Granulocyte colony stimulating factor. 

10. A system for individualized optimization of treatment of neutropenia, the 
10 system comprising: 

a system model comprising: 

a mathematical model of a neutrophil lineage; and 

a mathematical model of effects of treatment on said neutrophil 

lineage; 

15 a plurality of treatment protocols; 

means for adding individualized parameters to said system model, 

whereby said system model is modified based on said individualized 
parameters; and 

a selector for selecting an optimal treatment protocol from said 
20 plurality of treatment protocols based on modification of said system model. 
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higher efficacy can be obtained using 4-fold reduced total dose of TPO. 




TPO use in healthy donors: 
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Fig. 8: TPO given to healthy donors- Results of TPO clinical trials from recent research on healthy platelet donors, as 
25 compared to our computer simulation results. Arrows indicate the start of TPO treatment. (A) Comparison of experimental 
data from published articles 1 (black) and our model simulation (green), in both TPO was given as a single IV dose of 1.2 
Hg/kg on day 0. (B) Comparison of the same experimental data (black) and our proposed TPO administration protocol; the 
total dose in the simulated protocol was 0.3 jig/kg (blue). 



TPO use in patients receiving chemotherapy: 
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Fig. 9: TPO with chemotherapy- (A) Results of clinical trials from recent research on thrombocytopenia induced in patients 
50 receiving single carboplatin chemotherapy 2 on day 0 (black), as compared to our model simulation of these results (green). 
(B) The same experimental data (black); simulations of the same experiment, with addition of "conventional" TPO protocol 
of a single IV dose of 1.2 ug/kg on day 0 (olive); simulations of the same experiment under our proposed protocol that totals 
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APPENDIX C 

Optimizing cytotoxic drug delivery (administration/efficacy) for 
cancer patients 

Introduction 

Cancer is the second leading cause of mortality in the US, resulting in approximately 
550,000 deaths a year. There has been a significant overall rise in cancer cases in 
recent years, attributable to the aging of the population. Another contributing factor to 
the rise in the verifiable number of cases is the wide r u se of screening tests, such as 
mammography and elevated levels of prostate specific antigen (PSA) in the blood. 

Neither better detection nor the natural phenomenon of aging, however, can entirely 
explain the increase in new cases of tumors. Meanwhile, other cancers, like brain 
tumors and non-Hodgkin's lymphoma, are becoming more common. Their increase 
could refl ect changes in exposures to as yet unidentified carcino gens. Current trends 
suggest that cancer may overtake heart disease as the nation f s no. 1 killer in the 
foreseeable future. As gene therapy still faces significant hurdles before it becomes an 
established therapeutic strategy, present control of cancer depends entirely on 
chemotherapeutic methods. 

Chemotherapy is treatment withdrugs to destroy cancer cells. There are more than 50 
drugs that are now used to delay or stop the growth of cancer. More than a dozen 
cancers that formerly were fatal are now treatable, prolonging patients' lives with 
chemotherapy. 

Treatment is performed using agents that are widelyjion-ca^^ 
thatjiaye a high proliferation rate. Therefore, in addition to the malignant cells, most 
chemotherapeutic agents also c ause severe side-effec ts because of the damage 
inflicted on normal body cells. Many patients develop severe nausea and vomiting, 
become very tired, and lose their hair temporarily. Special drugs are given to alleviate 
some of these symptoms, particularly the nausea and vomiting. Chemotherapeutic 
drugs are usually given in combination with one another or in a particular sequence 
for a relatively short time. 

Chemotherapy is a problem involving many inter^l^^no^ 
operate on different organizational levels of the biological system. It usually involves 
genomic dynamics, namely, point mutations, gene amplification or other changes on 
the genomic level, which may result in increasing virulence of the neoplasia, or in the 
emergence of drug resistance. Chemotherapy may affect many events on the cellular 
level, such as cell-cycle arrest at different checkpoints, cell transition in and out of the 
proliferation cycle, etc. Chemotherapy may also interfere with the function of entire 
organs, most notably, with bone marrow blood production. In recent years molecular 
biology and genetics has made an important step forward in documenting many of 
these processes. Yet, for assessing the contribution of specific molecular elements to 
the great variety of disease profiles, experimental biology must be p rovided. wilLlPols 
that allow a ^ formal and systematic analysijs of the. intricate interaction between, the 
genomic, cellular and cell population processes in the host and in the disease agent. 
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This system is so complex that there is no intuitive way to know how small changes in 
the drug protocol will affect prognosis. But in spite of this intricacy, attempts to 
improve chemotherapy have been carried out by "trial and error" alone, with no 
formal theory underlying the application of specific drug schedules. Such an approach 
"is apt to result in no improvement, only discouragement and little useful information 
for future planning" (Skipper, 1986). 



1. General Description 

The product is a method for identifying cli nically feasible and efficacious optimal 
treatment protocols for anticancer drugs. We do this by using (1) optimization 
methods taken from Operation Research, and involving heuristic search algorithms. 
These are applied for (2) testing the performance of a very large number of drug 
protocols on (3) "virtual cancer patients," in the form of realistic computer 
simulations of detailed mathematical models of the involved dynamics. Our method 
is tKefirst which can identify optimal treatment protocols, without any 
constraint on the level of complexity of the mathematical models of the 
biological, clinical and pharmaceutical processes. This means that the 
description of the system can be as realistic as required. Previous optimization 
methods imposed constraints on the level of complexity of the mathematical methods, 
so that the simulated systems were always too simplistic to retrieve the clinical 
situation. 

The assumption that underlies this method is, that the effect of a certain drug on a 
certain cancer cells population depends on the overall population dynamics of the 
tumor. To describe these complicated dynamics we use a set of models thatjfollow 
up the changes in cancer cell numbers, spa tial distr ibution and other relevant 
phenotypic (e^eirieBc) factors, taking into consideration also the g enetic profile of 
th e patient , in cases where this one is available. To also allow for the use of cell- 
cycle-phase-specific drugs, we follow-up each drug susceptible body cell (normal and 
cancer) according to its cell-cycle status. Because any realistic description of such 
"virtual patient" is too complex to be attacked by analytical optimization method, then 
our introduction of the heuristi c optimization methods appro ach is a highly significant 
novelty. 

The project consists of 6 successive stages: 

1) Develo ping the appropriate m athematical models to describe the cancer cell 
populations with respect to the c ell-cycle as weH a^spatio-tempbfa l dyriaLmics. This 
model can be as complex as necessary for enabling an accurate computation of the 
number of cancer cells at any given moment. 

2) Inp utting individual parameters i nto the model and, hence, bridging between the 
model and accessible clinical data regarding profiles of the tumor and of the patient. 

3) Introducing drug characteristics (mode of action, effect and on tumor cells etc.), 
the pharmacolanefi the pharmacodynamics into the model, using experimental 
data about the relevant drugs. 

4) Combining this model _with addition al models that consider toxic effects of the 
drug on the host (e.g., granulocytic and thrombocyte lineages). 
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5) Defining a general fitness function f or given protocol, so that this system may 
serve as prediction tool for anticancer treatment (with an option of tuning by the 
user). 

6) Using sean^ijdg^ search through the possible 
protocols space in order to obtain the best/most suitable solution. 

The search process may be performed on several levels, starting from the most 
general (for a given type of disease), to the most detailed and patient-specific level 
(which will require relevant clinical data about the patient). In the case of a more 
generalized search, the certain region of protocols may be proposed as most suitable 
for given type/profile of cancer; in the case of specified search, this region can be 
narrowed according to additional patient data. 



2. Mathematical Model: 

2.1. Basic Tumor Layer 

The basic layer of the model incorporates a description of age distribution of cycling 
cells and number of resting (quiescent) cells. The term "age of the cell" here refers to 
chronologicalage startingjro m the conven tional beginningpoint of mitotic cycle. 

In what follows we give the definition of this basic model: 

- The whole cycle is divided into 4 compartments, or stages (Gi, S, G 2 and M). Each 
compartment is divided into equal subcompartments, where i subcompartment in 
each stage represents cells of age i in this stage (i.e. they have spent i time-steps in 
this stage). 

- Let {G/, 5, G2, M) denote the corresponding stages of the cell cycle, and Go - the 
quiescence. 

- Let 7jt denote the maximum duration of k** 1 stage in the cycle. 

- We denote as • t the time resolution of the model - discrete time steps. 

- x\ (t) is a function which represents the number of cells in stage k with age i into 
this stage, at time t to *+• Both time and age are measured in same units (hours). 

- Q(t) represents the number of resting cells at time t to *+• t. 

- Trans(k,i,t) represents the probability that a cell of age i in the stage k will move to 
the next (£+7) compartment. Cells entering the new stage always start from the first 
subcompartment, i.e. from i-l. This probability may change with time. 

By definition: \/k : Trans(k,T k ) = 1 That means that the cell cannot remain more than 
Tk in k compartment. 

- Let T 0 and T R denote the onset of Gl and the age at which the cell passes through R- 
point (restriction point) in Gl, respectively. The later is associated with commitment 
to complete the mitotic cycle. 

- The total number of proliferating cells P(t) can be calculated as follows: 

k=G\J y G2M\ /=1 J 



- In every time interval, quiescent cells may return to the proliferation pool. 
Alternatively proliferating cells may change their state to become quiescent if they 
are in Gl stage and at age /, where 7> i>0 . To describe this process we introduce the 
function Gy. o(ht) which describes the number of Gl cells in age i which become 
quiescent , during time interval [t, t+ • t]. This function may receive negative 
values, accounting for cells that return from resting to proliferation. 

As we assume that the exit to quiescence can occur only prior to the R-checkpoint 
(even in cancer cells), and that a resting cell which returns to proliferation enters the 
cycle at To, it can be stated: 

\/i>T R ,\/t:G^ 0 (i,t) = 0 

It must be noted, that this function is not dependent on i and t solely. Its value is 
determined according to current cell distribution and all the general parameters that 
characterize the described cells group. The same should be said about the values of 
Trans vector that can change during the history of given population. 

- The model traces the development of described group of cancer cells using given 
parameters, by calculating the number of cells in each and every subcompartment 
according to the following stepwise equations: 

jcf 1 (t - 1) • [1 - Transik , / - U - 1)] , \<i<T k 

T 

J [ jc/_, {t - 1) • Tmns{k ~lj,t- 1)1 / = 1 
jP 1 {t - 1) • [1 - TransiS, i-l,t- 1)], 1 < i < T Gl 

for k-GiM, k-l returns the previous stage (e.g. G2-1=S). 

2 [x J m it-D- TransiGl, j, t - 1)] • [1 - G^ 0 (i - 1, t - 1), i = 1 
>' it - 1) • [1 - TransiGl, i-lj- 1)] • [1 - G^ 0 (* - 1, t - 1)] , 1 < i < T Gl 

T 

l-Q^xiit-V-TransiMJJ-m, i = l 



(o = < 

x' it) = < 

x l it) = 

Gl 



- These equations make it possible to calculate the number of cells in each 
subcompartment at every time interval [t* t] starting from initial distribution (e.g. at 
time t=0). 

Remarks: 

1) As stated above, cells ages are measured in absolute time units, and therefore refer 
to the chronological age of the cell (as opposed to biological age, which refers to 
maturity of the cell and its units are relative to maturation rate that is different from 
cell to cell). As a consequence, in this model no cell can remain in the same age 
subcompartment after every time step. On the other hand, a fraction of the cells that 
leaves any subcompartment may be transferred to the first subcompartment of the 
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next stage, according to probability vector Trans(k,i,t). This fact enables us to account 
for variability of cycle lengths while retaining deterministic approach. 

2) The behavior of the cell population in this model is wholly controlled by two 
components: Trans vector, and GltoGO function. These two functions determine 
uniquely the outcome of every single time step, and, consequently the result over the 
long periods. We name them here "control functions". As we have mentioned, the 
values of these functions may be dependent not only on time and age of cells, but also 
on the current population state (or, generally, on the whole history of the population) 
as well as on the environment associated with given cells group. On the other hand, 
those parameters are similar for all the cells in the group. This implies that the model 
presented here is suitable for describing highly homogenous group of cells. As it will 
be claimed later, control functions values are determined by the parameters that 
describe the local conditions for given population. Therefore, we expect the basic 
layer of the model to give realistic description of uniform group of cancer cells for 
which environmental conditions and relevant biological properties are defined, in a 
way that will allow the construction of the control functions for this group. 

2.2. The General Tumor Model 

In the general approach the whole model is viewed as constructed from similar 
components, each of them derived from the basic structure described in the previous 
paragraph. Each component represents cells that are situated in the same 
environmental conditions and, therefore behaves similarly. We call them groups. The 
whole tumor is modeled as union of heterogeneous groups of cells, where the 
development of every group can be accurately predicted (when local conditions are 
known). 

This general model simulates progress of tumor in time by discrete steps. At each 
step the number of cells in each subcompartment of each group is calculated 
according to the previous state, parameters of tumor, drug concentration, etc. The 
parameters of the tumor must include all the information that is relevant to the 
prognosis. Some of these parameters are defined locally, e.g., those relating to the 
tumor 's geometry. For this reason the representation of the spatial structure will be 
included. 

The cells will be able to pass between the groups during the development of the 
tumor. This allows the representation of the changes in the local conditions, during 
the tumor evolution (e.g. forming of necrotic core, improvement in "living 
conditions" in vascularized regions, etc.). In addition, all the parameters of the tumor 
may change in accordance with the dynamics of the cancer. 

The calculation of the tumor development over time will be done by stepwise 
execution of the described simulation and can be used to predict the outcome of the 
treatment or in fitness function for search algorithms. 



2.3. From General to Individual Tumor Model 

When the general theoretical description of the model is accomplished, the model is 
fitted to represent the actual tumors. We render it patient-specific by adjusting all the 
parameters that determine the behavior of the modeled tumor to those of the real 
cancer in the patient's body. In order to accomplish this task we will establish the 
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connections between mathematical parameters (most of them will have direct 
biological implication) and every kind of data that is practically obtainable in the 
clinic. This individualization procedure includes three aspects: 

1) individual parameters of tumor dynamics 

2) individual parameters of specific drug pharmacokinetics for a specific patient. 

3) individual parameters of the dynamics of dose-limiting normal host tissues. 
These connections may be defined through research on statistic correlation between 

different parameters (including genotype-phenotype correlation), or using advanced 
biochemical research (which may establish quantitative relations). 

In the general model important dynamic parameters are estimated from experimental 
studies that are conducted in certain patient populations. Any of these parameters, 
when available on the per patient basis, can be individualized, while others left as 
population-based means. This approach allows continuous increase in the degree of 
individualization of the treatment protocols with progress in the technology. 

Defined so, the model will be able to give realistic predictions for treatment outcome 
either for specific patient or for a broad range profiles of patients and diseases. This 
tool can serve to perform the prognosis of an untreated cancer patient, or as a basis for 
treatment modeling. 

3. Introducing Pharmacology 

In order to simulate cancer treatment we add to this model the pharmacologic 
component. We model pharmacokinetics as well as pharmacodynamics for specific 
anticancer drugs. We begin by representing cell-cycle specific drugs, however the 
model implies no restriction on the type of drugs to be employed. 

We model the distribution of the drug in and around the tumor as well as in the 
blood (the drug kinetics). For this purpose, we use the suitable model, defining it 
precisely for every certain type of the drug. The concentrations of drug in the body are 
calculated at every time step in accordance with the drug administration specified by 
the protocol. 

The dynamics of the drug are represented through the direct influence of the drug 
on tumor cells. The effects on the proliferating cells are mostly blocking the cycle in 
different stages (which can be modeled as cell arrest) and cell death (immediate or 
after being in the block). Cell-cycle specific drugs have no direct influence on 
quiescent cells, but can affect them indirectly by killing proliferative cells and 
therefor changing local conditions. Where additional types of drugs added to the 
model, their effect on any kind of cells may be too modeled as killing certain fraction 
of cells (which is dose-dependent) or changing the behavior of the cells. 

Additional phenomena that may come out to be significant in drug kinetics and 
dynamics (e.g. rate of absorption by the cells, development of tumour resistance to 
specific drug, etc.) can be introduced into the model to make it as realistic as needed. 

The description of the drug in the model is done in terms of quantitative functions, 
which enable to calculate the drug amounts at certain locations and the tumor 
response to it at every time step. In the general case, these functions include 
parameters that depend on the specific data (drug type, body parameters, 
characteristic of the tumor, etc.) and can be determined in given situation. The relation 
between clinical data and these parameters can be established in the ways similar to 
those described for the cancer model. 

The combination of cancer model with the drug model described above makes it 
possible to predict the outcome of the treatment, given the relevant parameters for the 



drug, the cancer and the patient. Again, the prognosis may be made for specified cases 
as well as for broad profiles of patients or disease. This simulation also serves to build 
the fitness function used in search algorithms in optimization system. 

4. Combining with minimizing host toxicity. 

Although being accurate predictive tool, the model that represents chemotherapy of 
tumor alone cannot be used in optimization, for it posts no constraint on the choosing 
of the treatment. Actually, this model implies using as much drug as it is possible 
until the final elimination of the tumor; while in the live system the toxicity of the 
drug is the most important constraint limiting the treatment. In most cases of 
anticancer chemotherapy the dose-limiting toxicity is bone marrow suppression, the 
two most sensitive bone marrow lineages being granulopoiesis and thrombopoiesis. 
Accordingly, those two were chosen as an example and are modeled separately and in 
the similar way to predict the negative effect of the chemotherapy on them. These 
models reconstruct the damage caused by the chemotherapy to the bone marrow cells 
and the recovery of these lineages (also treated by specific growth factors). 

Thus, the whole system is able to predict the result of chemotherapy treatment for 
the tumor as well as for bone marrow cells, allowing the use of the protocols that 
combine anticancer drugs and growth factors for healthy cells. 

Chemotherapy toxicity to any other normal host cell population can be similarly 
taken into account, if it is proven to be relevant for dose and schedule optimization. 



5. Defining the Fitness Function 

Fitness function is an important tool in operation research. In the method developed 
by us it allows the comparison between a number of different protocols each one of 
them scoring differently with respect to various objectives that can be set by the 
developers or by the users and identifying the protocol for which we predict the best 
weighted score. Fitness function calculates for any given protocol its relative efficacy, 
thus enabling a definite decision of the best protocol from a given set of protocols. 

In different cases, different objectives can be formulated. There are 3 main settings 
in which such a model can be used: 
One) clinical practice- where objectives can change depending on type of disease, 

condition of the patient, purpose of treatment, etc; 
Two) pharmaceutical company- where objectives can be aimed at finding the 

therapeutic window and an optimal schedule; 
Three) scientific setting- where different and sometimes esoteric objectives 

can be aimed at. 

For any such case, certain fitness function can be formulated, reflecting all given 
requirements. Thus, in any particular case one can compare between different 
protocols and obtain the most suitable to his/her special purposes and needs. 

Using the fitness function we can: a) estimate the efficacy of a given protocol, b) 
search for the solution of an optimization problem, that is predict which protocol will 
be best of many potential protocols considered for curing/relieving the patient. 



6. Solving the Optimization Problem 

The optimization problem is stated using the described models: to find the protocol 
for drug administration (with option to growth factor administration) which 
maximizes the given fitness function. 

As explained above, we define the fitness function according to the user 
requirements. For example, the goal of the treatment may be defined as minimizing 
the number of cancer cells at the end of the treatment, minimizing the damage to the 
BM cells throughout the treatment or at its end, and curing the patient (where cure is 
defined precisely) as quickly as possible. Note that the fitness function may also 
include goals such as maximizing life expectancy, minimizing cost of treatment, etc. 

Generally, the aim of optimization is to find the best protocol, i.e. the protocol that 
generates the best value of fitness. Customarily, this is achieved by mathematical 
analysis. However, mathematical analysis is restricted to relatively simple models, 
whereas our models are very complex and, therefore, cannot be solved analytically. 
On the other hand, the practical purpose of the treatment is not to find the best 
possible protocol (i.e., the global optimum) but "only" one that will suit the user's 
objectives, even if its fitness is not absolutely the best (i.e., the local optimum). For 
this reason we can be satisfied with the solution that can be shown to promise that the 
patient will be cured. 

Hence, the optimization problem may be reformulated as follows: for given starting 
conditions, find the treatment protocol that will fulfill the user's requirements (e.g. 
curing a patient according to given definitions of cure) and subjected to given 
limitations (e.g. treatment duration, drug amounts, etc.). To this end it is unnecessary 
to find the global solution. It is enough to perform search, using search algorithms, in 
certain regions of the protocols' space, and find the local maxima of the fitness 
function. After determining these locally best protocols, we can verify that they serve 
one's objectives and check them numerically for stability. 

Such a strategy can be used for identifying patient-specific treatment, as well as in 
the general case, where only the profile of the disease and the drug specified. If more 
patient-specific data are supplied, the solution will be tailored more specifically. On 
the other hand, the optimization program could propose general recommendations for 
the protocol types for certain kinds of tumor. 



